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I. INTRODUCTION 



A. DIGITAL SIGNAL PROCESSING 

Digital signal processing has its roots in 16th century 
mathematics especially in the fields of astronomy arid com- 
pilation of mathematical tables. Today it has become a 
powerful tool in a multitude of diverse fields of science 
and technology. The application of digital signal process- 
ing varies from low-frequency spectrum seismology through 
spectral analysis of speech and sonar into the video spectrum 
of radar systems [l]. 

Rabiner and Gold have interrelated in their excellent 
book on "Theory and Application of Digital Signal Processing" 
(see Ref [2]) digital processing theory, 

with a variety of applications ranging from radar, sonar, 
communication, music, seismic and medical signal processing, 
and, with digital component technology. This technol- 
ogy is the main driving force for progress in this field 
as well as the general area of computer design. It is also 
observed by them that although the formulation of engineer- 
ing problems is often as vague as those of the "softer" 
sciences (such as anthropology, psychology, etc.) the execu- 
tion of these problems appears to depend on greater and 
greater accuracy and reproducibility. The capability of 
dig5_tal systems to achieve a guaranteed accuracy and essen- 
tially perfect reproducibility is very appealing to engineers. 
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A digital filter is defined in [3] to be a computational 
process or algorithm by which a digital signal or equiva- 
lently a sequence of numbers acting as input , is transformed 
into a second sequence of numbers or output digital signal, 
where, the term 'digital' implies that both time (the inde- 
pendent variable) and amplitude are quantized. A digital 
signal is said to be a discrete signal whenever the ampli- 
tude quantization step of a time sampled signal is zero. 

The field of one-dimensional digital filtering, which is 
outlined in Figure 1.1, encompasses recursive, nonrecursive 
and Fast-Fourier Transform processing. The terms recursive 
and nonrecursive, instead of HR (Infinite Impulse Response) 
and FIR (Finite Impulse Response) are preferred in the con- 
text of this realization oriented thesis. Kaiser showed in 
[4] that recursive processing is much more efficient than 
nonrecursive processing. Stockham's [5] method to perform 
fast convolution, which later became known as FFT method, 
improved the efficiency of nonrecursive techniques, so that 
comparisons in one-dimension are no longer strongly biased 
toward recursive methods [2]. 

One-dimensional recursive digital filter design depends 
stongly on the effective and well developed continuous fil- 
ter design theory. Moreover, the stability analysis of 
higher order one-dimensional recursive realizations can be 
solved by investigation of the root distribution of its fac- 
tored form representation. It is important to realize that 
continuous domain design techniques and the factorization 
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property (the fundamental theorem of algebra) exist only in 
one-dimension . 

Areas of applications of one-dimensional digital process- 
ing are listed in Figure 1.1. There are important military, 
and specifically naval applications in addition to those 
listed, such as missile and torpedo guidance, launcher con- 
trol system, fire control, combat information collections, 
dissemination and display. 

Digital filtering in several dimensions, which is over- 
viewed in Figure 1.2, has gained considerable importance, 
especially for the two-dimensional case. Figure 1.3 displays 
an example of two-dimensional low-pass and high-pass filter- 
ing. Until 1966, two-dimensional digital filtering was imple- 
mented by nonrecursive or convolutional techniques. In this 
method the output is the weighted sum of unit sample responses 
of all past input values. The serious disadvantage of the 
convolutional method is the requirement of a very large num- 
ber of arithmetic operations. The development of the Fast- 
Fourier-Transform (FFT) algorithm in 1966, reduced the number 
of arithmetic operations considerably and is used extensively 
today. Filtering via FFT is accomplished by computing the 
transform of the input function, multiplying by the filter 
transfer function, and inverse transforming the result. The 
recursive algorithm, which has in general an infinite impulse 
response [3], constitutes another technique for the realiza- 
tion of two-dimensional digital filters. 

Hall [6] compares the amount of computation required for 
the three filtering techniques. He shows that, in general, 
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FIG 1.3a EXAMPLE OF TWO-DIMENSIONAL DIGITAL FILTERING > [2] 
(ORIGINAL PHOTOGRAPH) 
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FIG 1.3b EXAMPLE OF TWO-DIMENSIONAL DIGITAL FILTERING, [2] 
(LOW-PASS FILTERED VERSION OF ORIGINAL PHOTOGRAPH) 
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FIG 1.3c EXAMPLE OF TWO-DIMENSIONAL DIGITAL FILTERING > [2] 
(HIGH-PASS FILTERED VERSION OF THE ORIGINAL 
PHOTOGRAPH) 



20 



the FFT and recursive algorithm are preferable to the non- 
recursive one in number of computations and storage require- 
ments. Hall demonstrates that the recursive filter algorithm 
constitutes the best method for large amounts of data, i.e., 
is the fastest and the cheapest. 

When processing two-dimensional data by recursive filters 
a fundamental problem exists due to the inherent feedback, 
namely, the problem of numerical stability. Since, in 
general, two variable polynomials cannot be factored into a 
product of real first and second order real coefficient poly- 
nomials, it is difficult to solve the stability problem for 
two-dimensional recursive filtering. Consequently, the 
majority of papers published in two-dimensional digital fil- 
tering deal with the design of nonrecursive filters which 
are inherently stable, i.e., [7] - [11]. 

There are several reports discussing two-dimensional 
recursive digital filters [ 12 ] — [ 1 5 ] . For example [ 1 2 ] and 
[15] formulate two-dimensional recursive filters by the Z- 
transform and linear difference equations, and although they 
investigate problem areas related to stability and realiza- 
tion, the majority of problems remain to be solved. 

The most important -applications of two-dimensional digi- 
tal filtering of digital data are in picture processing and 
geophysical data analysis. Picture processing can be cate- 
gorized into: 

1. Digital image restoration and enhancement, and 

2. Computer pictorial pattern recognition. 
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The methods of image restoration and enhancement are 
applied to invert degradations, such as aberrations, atmos- 
pheric effects, scanning, motion, and to manipulate images 
to improve viewing phenomena experienced by the human eye. 
Image restoration 'and enhancement has been used with great 
success in biomedicine [16], i.e., in extraction of quanti- 
tative information from X-ray films, chromosome counting, 
measuring the extent of atherosclerosis from arteriograms, 
in forensic sciences application, i.e., fingerprint image 
enhancement for automatic classification, and in astronomy 
[17], i.e., removal of turbulence from astronomical photo- 
graphy [ 1 8 ] . Pictorial pattern recognition by digital com- 
puter has gained great importance in satellite surveillance 
and military reconaissance . The application of two-dimen- 
sional digital filter permit the separation of different 
horizontal scans on magnetic and gravity maps in geophysics 
and can be used equally well for structural and topographic 
maps or for any other type of data which is available in 
the format of a planar grid [19], [4]. 

Although digital filtering in several dimensions is gain- 
ing importance in medicine, i.e. , heart volume measurements, 
in fire control problems, and to analyze complex electronic 
circuitry, there exists no general theory concerning struc- 
tures, analysis and design. 

These considerations are summarized as follows: since 
one-dimensional recursive digital analysis and design tech- 
niques depend strongly on the fundamental theorem of algebra 
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and on extensive utilization of continuous design theory, it 
represents a special, i.e., a non-generalizable field in the 
area of N-dimensional recursive digital filtering. Multi- 
dimensional digital filtering is well developed in two-dimen- 
sions but is based, due to the absence of recursive theory 
analysis and design tools, predominantly on nonrecursive 
filtering schemes. 

B. AREA OF INVESTIGATION 

The design of N-dimensional recursive digital filters is 
achieved by choosing filter coefficients to approximate given 
specifications and to evaluate the stability of the resulting 
filter. Several methods, which are discussed in Chapter II, 
are available to determine whether a multi-dimensional recur- 
sive digital filter is stable. They are based on a numerical 
root distribution investigation of the N-dimensional charac- 
teristic equation of H(z) and determine whether or not a fil- 
ter is stable.- There is no qualitative information extracted 
about how stable or unstable a system is or what coefficient 
changes of H(z) are necessary to make an unstable system 
stable. Moreover, these methods are complicated and result 
in time-consuming computer programs, which must be applied 
to each individual filter design variation. Due to these 
reasons, the design of multi-dimensional, recursive digital 
filters has not been achieved [2]. 

In the following chapters, a different approach to anal- 
ysis and design of N-dimensional recursive filters is pre- 
sented which is based on the Taylor expansion of the filter 

*The notation is discussed in Appendix A. 
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H(z). It can be shown that an N-dimensional recursive filter 
is stable if an only if the sum of all unit sample response 
entries converge absolutely, i.e., 




h 

( 0<n^<°°) 



where the h(n) are the coefficients of the Taylor expansion 
of N-dimensional transfer function H(z). 

Thus, conditions for stability of an N-dimensional recur- 
sive system HCz) can be obtained by observing the absolute 
convergency behavior of its Taylor expansion coefficients. 
Although this represents a rather complicated method to find 
stability conditions, it can be used to identify important 
unit sample response characteristics of certain two-dimen- 
sional recursive filters, which in turn can be generalized 
to N-dimensions . 

The relationship of the N-dimensional recursive filter 
coefficients to the unit sample response coefficients is 
investigated and methods to perform the approximation step 
of time-domain design are studied. Also, structural realiza- 
tion of N-dimensional transfer functions are searched. 

In summary a new design and implementation technique for 
N-dimensional recursive filters is developed. 

C. PREVIEW OF RESULTS 

The results can be grouped into techniques for stability 
analysis of N-dimensional recursive filters (Chapters V to 
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VII), three techniques to solve the approximation part of 
N-dimensional recursive filter design (Chapters VIII and IX) 
and, combining both previous results, design techniques are 
derived (Chapter X) . In addition, the implementation of 
N-dimensional recursive filters is considered and some new 
forms presented. These results are obtained using a compact 
vector index notation and a general theory defining charac- 
teristics of N-dimensional systems, as well as the N- 
dimensional Z-transform and some of its properties (see 
Appendix A) . 

The computations of N-dimensional Taylor expansion coef- 
ficients h(n) are discussed in Chapter IV. The standard 
derivative operator method requires complicated computations, 
and the convolutional or recursive technique requires the 
knowledge of all previous values to compute h(n) . Methods 
of combinatorial algebra are applied to derive two efficient 
nonrecursive formulas for Taylor expansion coefficients, of 
which one has high speed characteristics, i.e., requires very 
few multiplications and additions. The results of Chapter IV 
are used in Chapter V to solve the stability problem of the 
second-order (coupled), third- and fourth-order (uncoupled) 
two-dimensional recursive digital filter using Ultraspherical 
and derivative operator theory. 

In Chapters VI and VII, observations made in two- 
dimensional stability analysis of the previous chapter are 
generalized to cases of any order in two-and also in 
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N-dimensions j resulting in new sets of sufficient and 
necessary stability conditions. A general method to. derive 
necessary and sufficient stability conditions for N-dimensional 
recursive digital filters is presented and new stability 
conditions for the first degree and positive coefficient case 
in N-dimensions are listed. 

The important advantage of the new approach to the 
stability problem over the previously outlined numerical ones 
becomes evident if one considers that for a given type of 
filter transfer function the methods of Chapters V to VII are 
applied only once to obtain stability conditions in terms of 
coefficients of the characteristic equation of H(z) . These 
conditions can be used by the design engineer in a simple 
manner to test whether a specific filter is stable. They 
allow the identification of coefficients causing instability 
and can easily be interpreted to what extent coefficients 
must be changed to stabilize an unstable system. 

The nonrecursive combinatorial formula of Chapter IV is 
used in Chapter VIII to formulate rules for the propagation 
of the characteristic equation coefficients in H(z) in the 
corresponding unit sample response. Conversely, computer 
algorithms for the extraction of an exact all-pole or 
rational filter transfer function or a parallel arrangement 
of low-order all-pole or rational filter from a given unit 
sample response are presented in Chapter IX. All of these 
methods are new. The extraction of a parallel arrangement 
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of low-order sections can be viewed as an approximate method 
for partial fraction expansion in N-variables. 

In Chapter X, a time-domain design technique is presented 
using an extraction technique to approximate a given specifi- 
cation by a transfer function. The stability of the extracted 
transfer function is analyzed using the methods presented in 
Chapters V to VII. A measure of quality for the N-dimensional 
approximation techniques is introduced which compares the 
spectra of the given unit sample response with the spectra 
of the extracted transfer function. Finally, the implications 
of the N-variable partial fraction expansion technique on 
frequency domain design are discussed and areas of future' 
research are outlined. 

In summary a series method to analyze N-dimensional 
recursive digital filters is presented and new results 
regarding the stability of N-dimensional recursive digital 
filters are achieved. The series approach is used to formu- 
late three new transfer function extraction techniques which, 
combined with the stability and structural analysis, present 
a new time-domain design technique for N-dimensional 
recursive filters. 



«r. 
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II. METHODS TO DETERMINE STABILITY OF 
MULTI-DIMENSIONAL SYSTEMS 



A. INTRODUCTION 

Steiglitz [2l] has demonstrated by the construction of a 
specific isomorphism between the Hilbert Space of a complex 
valued Lebesque measurable function f(t) and the Hilbert 
Space of double-ended sequences of complex numbers { f (n )^__ OT > 
that the theories of one-diemnsional signal processing with 
linear time-invariant realizable filters are identical in 
the continuous and discrete case. It is interesting to note 
that optimization problems involving quadratic cost functions 
are also equivalent in the continuous and discrete time do- 
mains . 

The extension of this theory to several dimensions, 
though mathematically feasible, has little importance since 
a powerful continuous filter theory exists only in one di- 
mension. 

Another important property of one-dimensional theory is 
deduced directly from the fundamental theorem of algebra 
which states that every higher-order polynomial can be fac- 
tored into a product of lower-order polynomials. The impor- 
tance of the one-variable factorization property cannot be 
minimized. It permits the designer to investigate stability 
of a higher-order system by determining the factored form 
representation and applying simple stability criteria to 
each factor. If each and every factor is stable, then the 
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higher-order system is stable [l], [22]. The factorization 
property also constitutes a method of realization of higher- 
order filters by paralleling, or cascading lower order sec- 
tions, thereby reducing coefficient sensititity problems 
inherent in realizing higher order systems. Finally, an un- 
stable one-variable filter can be stabilized by cascading 
an appropriate all-pass filter [2]. This does not work in 
the analog domain because exact pole-zero cancellation is 
impossible. This technique would seem to be feasible in the 
discrete time domain except that roundoff effects may cause 
difficulties. It is shown in Appendix B that multivariable 
polynomials cannot be factored in general. The above are 
some of the major difficulties opposing the generalization 
of one-dimensional system theory to multi-dimensions. Some 
important conceptual ideas of one-dimensional analysis that 
can be generalized to several dimensions are presented in 
Appendix A. These include the definition of N-dimensional 
linear, shift-invariant, realizable, and separable systems, 
the N-dimensional Z-Transform and its properties , as well 
as the notion of the N-dimensional recursive equations and 
transfer function. 

The following paragraphs deal with existing methods to 
determine stability of N-dimensional recursive systems, H(z) 
(notation defined in Appendix A) which are obtained by apply- 
ing the N-dimensional Z-Transform to a linear, shift-invari- 
ant, causal N-dimensional difference equation. 
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B. METHODS TO DETERMINE STABILITY OF MULTI-VARIABLE SYSTEMS 



Shanks, Treital and Justice recently pulished a tech- 
nique to determine the stability of a two-dimensional recur- 
sive system, H(z 1} z 2 ), by which the unit disk in the z\- 
plane must be mapped into the z 2 -plane by solving the impli- 
cit two-variable denominator polynomial of H(z 15 z 2 ). A nec- 
essary and sufficient condition for stability is determined 
if the image of the z x -plane unit disk in the z 2 -plane does 
not overlap the unit disk in the z 2 -plane [ 12 ] . This method 
requires an infinite number of mappings and cannot, there- 
fore, be applied exactly. Huang [13] simplified the stabil- 
ity problem by the observation that 

A(z ! ,z 2 ) 

H( Zj, Z 2 = B ( 2i ,z 2 ) ’ 

where A and B are two- variable polynomials, is stable if and 
only if, 

1. The map of Sdi ={z i : | z i | =1} in the za-plane accord- 
ing to the equation B(z 15 z 2 ) = 0 lies inside d 2 = {z 2 :|z 2 |>l}. 

2. No point in d x = (z^lzj >1} maps into the point 
z 2 = 0 by the relation B(z 1 ,z 2 ) = 0. 

Although this procedure is finite, it requires the applica- 
tion of two bilinear transformations to formulate the prob- 
lem as discussed by Ansell [23] whose main contribution is 
to couple the use of a Hermite test for checking stability 
with a series of Sturm tests checking positivity. Anderson 
and Jury [24] modified Huang's method by outlining a pro- 
cedure which replaces the bilinear transformations by the 
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construction of a Schur-Cohn matrix and the checking for 
positivity of a set of self-inversive polynomials. It also 
replaces the Hermite test [25] by a Schur-Cohn matrix test 
and requires a series of Sturm tests or equivalently, a 
system of tests establishing the root distribution of a 
polynomial. These modifications represent a substantial 
reduction in computations as compared to Huang's method. 

Siljak [26] comments on [24] by proposing a substantial 
simplification of the number of Sturm tests which he states 
is desirable because a computer processing of the Sturm test 
is not available. However, a suitable computer algorithm 
for the Sturm test has recently been published by Singh and 
Panwar, [27] . The determination of stability of multivari- 
able filters is solved by Bose and Kamat [28] using a method 
outlined in [29] where the N-variable characteristic equa- 
tion of H("z) is reduced to several single- variable polynomi- 
als by an algorithm composed of a finite number of rational 
operations, followed by repeated factorization of the set 
of one-variable polynomials and back substitutions. This 
procedure is very complicated and the implementation on a 
general purpose digital computer may not be feasible [28], 
as recognized by the authors themselves. A slightly simpler 
method for determining the stability of two-dimensional 
filters has been offered very recently by Siljak C 3 0 ] in an 
algorithmic approach requiring the equivalence of two Marden 
tests [29] and checking positive definiteness with methods 
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developed in [3l], [32], and [33]. This last approach has 
been implemented by Schaldach [34] on a general purpose com- 
puter. Finally, Jury [36] presented an excellent summary 
on methods to determine stability of multivariable filter 
using the inners approach. 

C. SUMMARY 

The above outlined methods to determine the stability of 
multivariable characteristic equation of H(z) to identify 
roots not satisfying 

n 

Pi |z. | < 1 (2.1) 

i=l 

where n iz.i < i , (| Zl i < i)n (i z 2 1 < i)n...n(iz n i< i). 

They also have in common the property that they are highly 
theoretical and very difficult to apply [2] and of limited 
practical value to the design engineer. It is necessary for 
analysis and design purposes of multivariable systems to 
identify not only whether or not a system is stable, i.e., 
that all roots of the N-variable characteristic equation 
satisfy eq (2.1), but also what coefficients of an unstable 
transfer function cause the instability and how the coef- 
ficients must be changed to achieve stability. This latter 
question can only be solved definitely utilizing previously 
considered methods for the bilinear two-dimensional case, 
for example: 
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H(Z 1} 2 2 ) 



1 - a in z - 1 - ~ - - 1 - ~ - - 1 - - 1 



'10 1 



~ a 01 Z 2 



‘ «11 Z 1 



( 2 . 2 ) 

for which stability conditions are derived in [12 ], [13], 
and [24]. This result will be discussed further in Chapters 
V and VII. 
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III. STRUCTURES AND STATE SPACE FORMULATION 



FOR N-DIMENS TONAL DIGITAL FILTER 
A. INTRODUCTION 

There are several ways of realizing a given digital filter 
transfer function. The choice of structure depends on finite 
wordlength effects in practical hardware realizations on 
complexity of computation. 

One-dimensional canonical structures are most economical 
and fastest, yet they may have undesirable finite register 
length effects. In this chapter a one-dimensional direct 
form filter realization is introduced designated by direct 
form 4 to distinguish it from other direct forms which 
corresponds directly to hardware implementation and by 
realizing the numerator and denominator polynomial separately 
in cascaded form. The direct form 4 is deductively used to 
derive direct form and canonical structures in N-dimensions , 
which have considerable importance, since the nonfactora- 
bility characteristic of multivariable polynomials eliminates 
the possibility of realizing high-order multi-dimensional 
filter by paralleling or cascading low-order standardized 
sections, except for the separable case. 

The state space .approach provides a means for character- 
izing complex systems in a simple manner. The complexity of 
manipulation is decreased and the easy implementation on a 
general purpose digital computer minimizes computational 
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difficulties. The increase of the number of state variables * 
the number of inputs, or the number of outputs does not in- 
crease the complexity of the equations. In fact, the analysis 
of complicated multiple- input-multiple-output systems can be 
carried out by procedures that are only slightly more compli- 
cated than those required for the analysis of systems of first- 
order equations. 

The one-diinensional state variable approach to digital 
filtering was presented by Grindon [37] in 1970 . Since then it 
is commonly used in one-dimensional digital filter analysis. 

For example, Parker and Hess [ 33 ] used the state space approach 
to derive twenty four second-order canonic digital filter sec- 
tions . 

Parker, Girard and Souchon [ 42 ] used state space methods 
to express the co-variance matrix of state errors and the 
variance of the output error in terms of the correlation ma- 
trices of the quantization error sources. This method allowed 
the evaluation of correlation effects between error sources in 
the comparison of different filter structures. 

Parker and Yakowitz [4 4j used the state space approach to 
compute quantization error bounds. 

It can be summarized that the state space approach is a 
powerful and commonly used method to analyze one-dimensional 
digital systems. Since the following chapters are mainly con- 
cerned with N-dimensional filter theory, an extension of the 
one-dimensional state space analysis to N-dimensions is devel- 
oped in Section C of this chapter. 
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B. STRUCTURES OF N-DIMENSIONAL DIGITAL FILTER 



1. One-Dimensional Digital Filter Realization 

The rational transfer function (A. 10) has for the 
one-dimensional case, i.e., N = 1, the form: 



H( Zj ) 



L 



^ b z 
/ . Zl 



£= 0 

M 




m=l 



Y(Zj) 

X ( z t ) 



(3.1) 



The time domain representation is obtained by applying the 
inverse Z-Transform to eq (3.1) which results in: 



y(n, ) 



M 

Z a y(n.-m) 
m J i 

m=l 



•i 

o - n 



3^ x(n x -£.) 



(3.2) 



The realization implementing eq (3.2) directly, is known as 
direct form 1 and is shown in Figure 3.1. 

Equation (3.1) can be rewritten in a slightly different 
form by introducing a new variable, W(z x ), such that 



W(zi) Y( z ! ) 

H(Zj) = X(Zj) * W ( z j ) 



(3.3) 



The corresponding set of difference equations consists of: 



i 
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(3.4) 



w(n ! ) 



x(n ! ) 



M 




m=l 



and 



Y (n ! ) 



L 

^ 6 ^ w(n !-£) 

s,=o 



(3.5) 



Equations ( 3 . 4 ) and (3.5) are realized in the direct form 
2 as shown in Figure 3.2 which can be visualized as a cascade 
realization of two digital filters, realizing the denominator 
and numerator polynomials, respectively. Figure 3.2 can be 
redrawn as shown in Figure 3.3. The resulting realization 
is called direct form 3 or canonic form, since it represents 
a structure with the number of delays equal to the order. 

At this point, an additional realization will be intro- 
duced, the direct form 4, which combines the characteristics 
of the direct forms 1 and 2 by having only two entries in 
each summer, which corresponds directly to hardware implemen- 
tation and by realizing the numerator and denominator poly- 
nomial separately in cascaded form. The usefulness of the 
direct form 4, which is shown in Figure 3.4 and which re- 
places the direct forms 1 and 2 , will become evident in the 
following section. 

There are several other methods to realize a one-dimen- 
sional order digital filter. For example, the cascade 

form of first and second order sections H(Zj), as shown in 
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Figure 3.6, where 



H(Zj) 



k 

n 

i = l 



Hi (Zj) 




and [x] = largest integr 



(3.6) 
< x 



Each component H^(Zj ) can be realized in one of the above 
outlined forms. 

Partial fraction expansion methods applied to eq (3.1) 
lead to a parallel from realization of first and second 
order sections (see Fig 3.5) i.e., 



H( z j ) = C 



I 



i = l 



Hi ( z ! ) 






(3.7) 



Other forms of realization include hybrid structures, i.e., 
parallel-cascade forms, transpose configurations, which can 
be obtained for all previously outlined structures by revers- 
ing the direction of signal flow and by interchanging all 
branch and summing modes, wave-digital filter [46] to [50] 
and continued fraction expansion filter [51], 

The direct forms discussed above and the series/parallel 
arrangement of lower order sections are by far the most often 
used ones in computer simulation and in digital hardware. 



2 . Deductive Approach to Develop N-Dimensional Structures 

The fundamental theorem of algebra does not extend 
to the multi-variable case which is proven in Appendix B. 
Consequently, the series and parallel realization schemes 
cannot be generalized to N-dimensions . The following 
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paragraphs develop a generalization for direct forms to two- 
dimension and then to N-dimensions . 



The numerator of the two-dimensional transfer function 



Li L 2 



I I V 



-£, -£,. 

£ 2 Z > lz 2 



H( z 1 , z 2 ) = 



£ 2 = 0 £ 2 = 0 



Mj M 2 



( 3 . 8 ) 



z -m, -m 2 

m 2 1 2 



“I I-. 

m x =0 m 2 =0 
(m 1 +m 2 ^ 0 ) 

can be rewritten as one-variable polynomial in z 2 , weighted 
by coefficients which are functions in z 2 , i.e. , 



N(z 2 ,z 2 



L 2 / L i 

'■Z{l 

£ 2 - 0 ■ ' £ i “0 



£ i \ _ Ho 



6 £ 2 £ 2 Z i M 2 



( 3 . 9 ) 



4 / 

Ik 

£,=0 ' 



. (Zl) 1 22'^ 



( 3 . 10 ) 



The same reasoning applied to the denominator polynomial 
leads to 



M, / M 2 



D(z ljZa ) = 1 - ]T ( 

m 2 =0 ' m 2 = 0 



m ^ \ — m 2 

*2 



( 3 . 11 ) 
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1 



(3.12) 



^2 . 

- Y ( i w zi) ) z2 ~ m2 

m 2 = 0 

where the sum of the indices is mi + m 2 ¥■ 0 . To obtain 
the direct form 3 representation, a new function, W(zj,z 2 ), 
is introduced: 



W( Zi , z 2 ) Y ( zi, z 2 ) 
H(z 15 z 2 ) - x(z 15 z 2 ) ' W ( Zi , z 2 ) 



(3.13) 



where each factor can be written using eq (3.10) as 



W(z 15 z 2 ) 



L 2 

■U\^) z 2 ^ 2 x(z i’ zz) 
£ 2 = 0 



and eq (3.12) as 



(3.14) 



Y( Zj , z 2 ) 



M 2 

W(Zj,zp + ^ ( D Mi m 2 z 2 ' m 2 Y(z 1? z 2 ) 

m 2 =0 



(3.15) 



The realization of eqs (3.14) and (3.15) is shown in Figure 
3.7 using compact notation, and in Figure 3.8 where 

an< ^ m are implemented for each (L : £ 2 ) and 

and (t^m 2 ) by one-dimensional direct form 4 realizations. 
The direct form 3 (canonic) will be derived in the next 
section. The direct form 4 realization of H(z x ,z 2 ) is 
shown in Figures 3.10 and 3.11, using compact notation and 
one-dimensional direct form 4 implementation of each N T 0 

1j j 
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and , respectively. The generalization to N-dimensions 

follows directly where the N-dimensional transfer function 
is written in the numerator and denominator as a one-vari- 
able summation which is weighted by functions of (N-l) vari- 
ables , i . e . , 



The direct forms 3 and 4 in N-dimensions are shown using com- 
pact notation in Figures 3.12 and 3.13, respectively. 

Example 3.1: 

The direct form 4 realization of a full quadratic trans- 
fer function in three-dimensions is shown in Figure '3.15. 
Structures of two-dimensional third-order and fourth-order 
uncoupled, as well as second-order and fourth-order coupled 
filters, which will be investigated in detail in the next 
chapters, are presented in Appendix G. 

A separable system is defined in Appendix A(4) by 



L 



'N 



H(zi, . . . ,z N ) 




(3.16) 



N 

h(n) = n h.(n.) 
i = l 1 1 



(3.17) 



which is shown in Appendix A(5) to imply: 
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N 

H(z) = n H . ( z . ) (3.18) 

i=l 1 1 

But H(z) can be realized by a cascade structure of one-dimen- 
sional transfer functions in each dimension, as shown in 
Figure 3.16, by techniques discussed in Paragraph 1 of this 
section. 



C. STATE EQUATION FORMULATION 

1. State equation formulation in two-dimensions 

The two-dimensional transfer function H(z 15 z 2 ) can be 
written in long hand notation as 



Li L 2 



I L\ 

£1=0 £2=0 



- 2-1 - 2-2 

i z Z I 



H(z,,z 2 ) = Mi 



-7 7“ z 

Z_, m i m 2 1 



- m i Z " m 2 



mj=0 m 2 = 0 

(m 1 +m 2 ^0 ) 



(3.19) 



If the function W(z x ,z 2 ) is introduced similarly, as in the 
previous section, then 



W(Zj,z 2 ) Y(Zj,z 2 ) 
H(Zj,z 2 ) - u( Z j J z 2 ) W ( Zj , z 2 ) 



Each factor can be written as: 
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L L 
1 2 



Y ( Zj , Z 2 ) = Y 

£i =0 ^ = 0 



^ Zj’ 5,1 Zj"^ 2 W(ZjpZ 2 ) (3.20) 



and 



Mj M 2 

W(Zj,z 2 ) = U ( Zj , z 2 ) + ^ Y t a m ^ m ^ z i~ mi z 2 ~ m2 W(z t ,z 2 ) 



mj=0 m 2 =0 
(mj+1112 iO ) 



(3.21) 



The inverse Z-transforms of eqs (3.20) and (3.21) are 



y(ni,n 2 ) 



Li 

= V 

=0 

x 




w(n 1 -£ 1 ,n 2 ~£ 2 ) 



(3.22) 



and 



w(n 15 n 2 ) 



M x M 2 

a w(ni-m 2 ,ni-m 2 ) 

mjmj 1 ’ 

nq^O m 2 = 0 
Cirij+m^O ) 



u(ni,n 2 ) + 



Z Z 



(2.23) 



Without loss of generality, l,£ is chosen equal to for all 
i. The structure realizing eqs (3.22) and (3.23) is shown 
in Figure 3.9. It is noted that the number of delays in 
Figure 3.9 equals the order of H(z p z 2 ). The structure is, 
therefore, by definition, canonic. 
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Definition 3 : 



The state Xp^(nj,n 2 ) is defined to be a number which 

S "t St 

enters at the time instant (n x ,n 2 ) the (p+1) , (q+1) delay 



m z, , z 2 



-1 direction, respectively. Following definition 
3 the states are marked in Figure 3.16. Comparing eqs (3.22) 
and (3.23) with the structure in Figure 3.16 easily identi- 
fies the following equalities: 



X 00 (n i’ n 2 ) 


= w(n l5 n 2 ) 


X 10 ^ n i ,n 2 ^ 

• 


= w ( nj - 1 , n 2 ) 


• 

X MjO (n i ,n 2 ^ 


= w(n 1 -Mj , n 2 ) 


X 01 ^ n i ,n 2 ^ 


= w(n x ,n 2 -l) 


^ll^ n i’ n 2 ^ 


= w(n!-l,n 2 -l) 


X M1 (n ,i’ n 2 ) 
1 


= w(n 1 -M 1 ,n 2 -: 


X 0M 2 ^ n i ,n 2 ) 


= w(n 1 ,n 2 -M 2 ) 


X lM 2 (n i’ n 2 


) = w(nj-l,n 2 - 


X M 1 M 2 (n 1 qn 2 ) 


= w(n 1 -M 1 ,n 2 -M 
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From Figure 3.16 the relationship between the state vari- 
ables is identified as: 



Xqo ( n p n 2 ) - a i0 


X 00 (ni_1,n 2> 


+ 


a 20 


X^C^-l,^) 


+ . . . 


+ a 30 


X 0 0 > n 2> 


+ 


a 01 


x 00 (n i’ n 2 _1) 


+ . . . 


+ a ll 


X Q1 (n 1 -l,n 2 ) 


+ 


a 21 


X 11 ^ n i _1 ,n 2 ^ 


+ . . . 



+ a 0M 2 X 0 (M 2 a ) + a iM x qm ^ n i~l» n 2^ + '* 



+ a 2M X 1M + • • • 

2 2 



X 



10 



(n 1 ,n 2 ) = Xq Q (n 1 -l,n 



) 



^01^^1 5 ^2^ X QQ ^^1 5 ^2 ^ 

X 11 (n 1 ,n 2 ) = X Q1 (n 1 -l,n 2 ) 



X 0M^ n i’ n 2^ = X o (M 2 -l) ^ n i ,n 2 
^1M l > n 2 ^ - XqM (n 1 -l,n 2 ) 



These equations are written in matrix form as follows: 



where 



— ( n i’ n 2 ) = 



X 00 X 10 x 20 *•' X (M r l)0 X 01 X 11 X 21 



X (Mj-l) 1 X 0M 2 X 1M 2 X (M 2 -1)M 2 | 



rt 









1 

l°l 



_ r i — t -i 

Aj | 0 |. . .| 0 



r 



R 



i i 
i i 

i°i 



i- 

r 



t- 



L 



£1 

0 



I 0_l 

I I 
I I 
l°l 
I I 



ro 



R 



I . I 
I • I 

_J •_] 

I I 
I 0 I 



r. 



• LrL’-iJ 

0 ^ 



10 



. . . ( 0 



...10 



. . ( 0 






I 0 



! l_l LLJ 






R 



where 



A* — I Ot - Ot r \ • C l rt • ... Otw . 

D L lo 2] 3] Mj]J , 



1 

0 

0 



0 

1 

0 



0 

0 

1 



0 0 0 



00 

00 

00 



10 



, which is (Mj-1) X Mj 



and 0 is the zero matrix. 
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\ 

Q( 0 ) 


1 

'. . . 


1 

1 


Q(0) 
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L 
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where 



Q(p) = 



P 


0 


0 ... 


0 


0 


0 


0 ... 


0 


0 


0 


0 ... 


0 


0 


0 


0 ... 


0 



which is X M l 



B = [l 00. ..0 000. ..0 ...000 ... 0] 



The output equation is similarly formulated in matrix form: 



y(n lS n 2 ) = C a X(n 1 -l,n 2 ) 



C 2 X(n l5 



n 2 ~1 ) 



Du(n l5 n 2 ) 



(3 



where X is defined as before, and 



• 25 ) 
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Bio + a io Boo t 


C t r 
_1 2 


801 + a oi ^°° 


B 2 o + a 20 Boo 




0 


B 30 + Boo 




0 


B 00 




^02 + a 02 ^00 


^21 + a 21 ^00 




0 


^31 + a 3 1 ^0 0 




0 


3 i(M 2 -1) + a KM 2 -l) 


* 


3 0M 2 + a 0Mi, 3 oo 


S2 (M 2 -1) + a2 (M 2 -l) 3 °° 




0 


B3 (M 2 -1) + a 3(M 2 -l) S 00 




0 


8 w + ol 8 

1M 2 !M 2 p oo 




0 


6i m ai M 2 3)0 




0 


3l M 2 ai M 2 600 




0 



and D = B 00 

2 . State equation in N-dimensions 

The procedure for obtaining state equations in two 
dimensions can quite readily be extended to N-dimensions, 
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x. e. 



5 



the equations for y(n) and w(n) are 



y(h) 



I 



Sj w(n-£) 



( 0 ) 



and 



w(n) 




w(n-m) 



(3.26) 



(3.27) 



The state and output equations 

X(n x , . . . ,n N ) = X(n x -1 ,n 2 , . . . ,n^) + A 2 X(n x ,n 2 -l ,n 3 . . .,n^ T ) 

+ ••• + — ^ n i 5 • • • > n N ~l) + Bu(n x , . . . ,n N ) (3-28) 

Y (n x , . . . ,n^) - C^ X( n i“l >n 2 , . • . , n^j ) t C 2 X(n,,n 2 — l,n 3 ,. . . , n ^ ) 

+ ... + C^j X(n x , . . . , n^-1 , Nj^— i ) Du (n j j • . . > n^ ) 

(3.29) 

D. DISCUSSION OF CANONIC SECTIONS IN N-DIMENSIONS 

It was stated in Section B of this chapter that a one- 
dimensional high-order transfer function can be implemented 
in general as a parallel, cascade or hybrid arrangement of 
second-order sections. For reasons of noise generation, 
coefficient sensitivity, economics, logistics and others, 
these low-order sections are very important. The canonic 
realizations, for example in Fig 3.3, which constitute 



the subclass of second-order systems having minimum number 
of components, are widely used today [2], [ 4 0 ] , and [ 4 1 ] . 

The non-factorability of multi-variable polynomials and 
the fact that no stability conditions except for the two- 
dimensional bilinear case (eq (2.2)) were available up to now 
explains why recursive digital filter design has not been 
achieved in two- or higher-dimensions. 

In" Chapter X, N-dimensional time domain design techniques 
are presented which result in high-order N-dimensional trans- 
fer functions or a parallel arrangement of low-order sections. 
The order of these parallel sections depends on stability 
and design considerations. It was shown that factorization 
or partial fraction expansion cannot in general be achieved. 
However, approximations (See Chapter IX) are possible, so 
that the consideration of canonic sections is of value. In 
Appendix F, a state space approach is proposed to derive 
multi-dimensional canonic structures for the general second- 
order two-dimensional case. The derivations of canonical 
sections in general are not further pursued, but are called 
to the attention of future researchers in the multi-dimen- 
sional digital filter area. 

E . SUMMARY 

In this chapter an N-dimensional structure for the imple- 
mentation of digital filters was presented. The direct form 
which combines features of the direct forms 1 and 2, was 
used to develop compact realization of N-dimensional filter. 
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The significance of these N-dimensional realization schemes 
stems from the fact that realization schemes can be used 
directly to implement hardware, where components of the 
realization are replaced by summers, multipliers, and shift 
registers . 
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FIG. 3.1 • DIRECT FORM 
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FIG. 3.2: DIRECT FORM 2 
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FIG. 3.4= DIRECT FORM 4 
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FIG. 3.7= DIRECT FORM 
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FIG.3.8: DIRECT FORM 
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FIG. 3.9 : DIRECT FORM 3, 2D 
(Canonic) 
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FIG.3.IO : DIRECT FORM 4 , 2D 
(Compact Notation) 
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FIG.3.II: DIRECT FORM 
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FIG.3.13: DIRECT FORM 4, ND 
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Fig. 3.14= BLOCK DIAGRAM REPRESENTATION OF A 
SEPARABLE SYSTEM. 
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(STATE ANALYSI 




IV. SERIES METHOD TO DETERMINE STABILITY OF 
AN N-DIMENSIONAL RECURSIVE FILTER 



A. INTRODUCTION 

The approximation part in the design of recursive filters 
is carried out by choosing filter coefficients to approximate 
a given specification in time or frequency domain and evalu- 
ating the stability of the resulting filter transfer func- 
tion. The generalization of one-dimensional design tech- 
niques for recursive digital filters, in particular stability 
analysis, to the multivariable case is not possible because 
the fundamental theorem of algebra which permits one to in- 
vestigate stability of higher-order systems by applying 
simple conditions to the roots of the factored form repre- 
sentation, cannot be extended to several variables. 

Time domain design techniques are developed 
in Chapters VIII, IX, and X. At this point, however, it is 
of interest to note that the stability investigation is a 
decisive part of the design procedure. From the design 
engineer's point of view, it is very impractical if the 
stability of each individual filter design variation must be 
determined using complicated methods which result in time 
consuming computer use, and when the filter coefficients 
causing instability cannot be identified and, consequently, 
no knowledge is obtained to what extent coefficients must be 
changed to stabilize an unstable system. Due to this fact, 
the design of multi-dimensional recursive filter has not 
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been achieved. In fact, there is not even a simple design 
method available [2]. 

In the following Chapters, an approach to the solution 
of the stability problem is presented which is not based on 
the N-dimensional roots of the characteristic equation. In- 
stead, it employs the N-dimensional Taylor series expansion 
of a transfer function whose coefficients are values of the 
unit sample response, h(n), i.e., the Taylor expansion of 
H(z) is: 



H(z) 



^ h(n) z n 
n 

( 0<n^<°°) 



(4.1) 



It is known that an N-dimensional system is stable if every 
bounded input produces a bounded output (= BIBO stability). 
For linear, shift-invariant recursive digital filter a nec- 
essary and sufficient condition for BIBO stability is 



S 




n 

( 0<n^<°°) 



< oo 



(4.2) 



Proof : See Appendix C. 



To identify stability conditions or conditions for abso- 
lute summability of the N- variable sequence of eq (4.2) the 
general term h(n) must be identifiable. Four formulas are 
presented in this chapter for the computation of h(n) which 
are based on derivative, recursive, nonrecursive, and non- 
recursive combinatorial techniques. Both nonrecursive 
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formulas express h(n) explicitly in terms of the N-variable 
characteristic equation coefficients. Several examples and 
a comparative evaluation of the formulas demonstrate the 
perspicuity of the nonrecursive combinatorial method. 

The above computational techniques are employed in' Chap- 
ter V to compute unit sample responses for the bilinear third 
and fourth-order (uncoupled) two-dimensional transfer func- 
tions. By application of derivative operator and Ultraspher- 
ical polynomial methods, necessary and sufficient as well as 
sufficient stability conditions are identified for the above 
cases. The concept of a "stability region in the coeffici- 
ent space" is introduced and pictures of the bilinear filter 
stability region are shown. 

In Chapters VI and VII, sufficient conditions and neces- 
sary conditions for the general N-dimensional filter, and 
methods to compute necessary and sufficient conditions for 
the general N-dimensional transfer function are derived. 
Several examples are presented to demonstrate the application 
of the above methods . 



B. TAYLOR EXPANSION FORMULA IN N- VARIABLES 

The Taylor series expansion formula for a function, 
F(z), of N-variables is given in [52] as: 




+ ... + 




F( z ) | (4.3) 

( 0 ) 



k ’=0 
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where 



( 



z, - — + 
1 3Zi 



J 2 3 z. 



J N 3z 



-Y 

N/ 




+ z 



N 




+ 



+ 



N 




(4.5) 



denotes a differential operator. Since the expansion will 
always be taken about the origin, we could also term it 
He Laurin expansion. Following [53], eq (4.3) will now be 
reformulated into a more intelligible form. 



The powers of any term, 



( z ‘ 3 0 '"'( Z 4 ;) N 



(4.6) 



of the expanded form of the product in eq (4.5) will satisfy 

ij + . . . + i^ = k. The number of times the term of eq (4.6) 

occurs, equals the number of different permutations of the 

* 

k factor. Thus, 




where i = [i 15 ..., i^] assume all positive values satisfying 
the relation ij . . . + i = k. 



* 




k! 



V 
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Equation (4.3) is rewritten, using eq (4.7), as:* 



F(z) 



.7>E(. * . ) f -A 

h? ’ V 1 !’ ' • • 5l N /3Z J . ..3z n 



XT F.(z) 1 z 
' N (0) 



l (4. &) 



(i + . . . + i^=k) 



which is equivalent to 



F(z) = 




js lj + . • • + i jq __ 

2-L — F(z) i Z 

azj 1 ... 3 z n 1n (0) 



-it it 

lj. ... J.JJ . 



(O^i^l 00 ) 



(4.9) 



Also, if we define 



.(i 1? . . . ,i N ) ( z) = 



9 i 1 +. . . +iN 
azi 11 ... z n x n 



F(z) (4.10) 



then eq (4.9) becomes 



F(z) = V) -* — - — =— - f (:Li ’ * ' ’ ’ 1n) (z) | z 1 (4.11) 

// . * * ‘ 1 N* (0) 

*r 

l 

( 0<i^<°°) 

This form of the Taylor expansion of a function F(z) will be 
used in the subsequent derivations in this chapter. Its 
application will be demonstrated in the following short 
example . 

Example 4.1 : 



The first few terms of the Taylor expansion of a third- 
order two-dimensional transfer function, 




1 



(4.12) 



L(z l5 z 2 ) = 



" a io Z 1 “ a 0 1 Z Z 1 



can be expressed, using eq (4.11), as: 



oo oo 



L(z„z 2 ) = 2 ^ 2^^-r^-r L <1 >’ 1 2 ) (z 1 ,z 2 ) I 



Z h z h 

Zj Z 2 



ii = 0 i 2 = 0 



( 0 , 0 ) 



= 1 + -t— L(z ,z ) | Zj + L(z ,z ) | z 2 + 

3z i (0,0) 9z 2 (0,0) 



+ sl - 37“ L( z , z ) | 2 z + L(z ,z 2 ) | z 2 z+... 

3z i 3z 2 1 2 rn. n) dz i 3z i rn.m 



3 Z 3 



( 0 , 0 ) 



( 0 , 0 ) 



1 + a 10 z, + a Q1 z 2 + 2d 10 a 01 z,z 2 + 



+ ( 3 a io 2 a oi + 2 a oi “20) z i 2 z 2 + 



(4.12a) 



C. GENERAL TERM OF THE N-DIMENSIONAL TAYLOR EXPANSION 

It is useful for later derivations to formulate the Tay- 
lor expansion of the N-variable transfer function H(z), i.e., 



H(z) = 



1 - Q (z -1 ) 



(4.13) 



1 



in terms of z . This will be possible if the expansion of 



L( z) = 



I - Q ( z ) 

can be found in terms of powers of z . 



( 4 . 15 ) 
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In three of the following sections it will be shown that the 
coefficients of the Taylor expansion of L(z) can be explic- 
itly derived in terms of the Taylor coefficients of Q(z). 

1. Derivative Formula 

If the Taylor expansion formula of H(z) is written 
as 




H(z) = )) h(i) z 



-l 



(4.16) 



( 0 < ij <po ) 



then by eq (4.11), h(i) is equal to 



h(i) = r— r H 



~ — 3 — r h^ 1 *' “ ,lN )(z) 



1 N ‘ (0) 



(4.17) 



2. Recursive Formula 



Theorem 4.1 : If the Taylor expansion of Q has the form 




i 

( 0 < ij <o° ) 




(4.18) 



and the Taylor expansion of 



L( z ) 



1 

1 - Q(Z) 



(4.19) 



is 
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(4.20) 



L( z ) 



= S] b (i) z " 



( 0 < j . <oo) 

“ 1 “ 



then 



’(I) 



= 2] a <i-J 



j) b (j) 



(4.21) 



D<i 



Proof : It follows directly from eq (4.15) that 

L(z) - L(z) Q ( z ) = 1. 

Substituting eqs (4.18) and (4.20), it follows that 



^ b d) z± 



At b (D 



a (D 2 = 1 



1 



When i = 0, since a^^ = 0 it follows b (Q) = 1* For all 
values of i ¥■ 0 it follows that 



lL b a> z ' 



b (D zI }) a (i) zl 



As a result 






3<i 



a (i-j) b (j) 



7 ^ 



Example 4.2 : 



The first few coefficients of the Taylor expansion for 
L( z i , z 2 ) , where 



L(z ,z } = 

1 " a io 



\ ~ a 0 1 Z 2 a 2 0 



(4.22) 



are computed by eq (4.21), where the only nonzero coeffici- 
ents of Q(z 15 z 2 ), which is defined by eq (4.18), are: 

a, , a ^ = “io ; 



L (i,°) 



(o,0 



— ot 0 1 , a (2 0 ) _CX20 



b , . . \ 

(i i , I2 ) 



x 1 1 2 



V V 

j 1 = 0 j 2 = 0 

J 1+0*2 ¥■ ii+i2 



It follows : 



r — 1 



( 1 , 0 ) / / a (l-jj,0-j 2 ) 

j-o j 2 =o 
3a + 3 2 ?t 1 



1 = a 



10 



(0 , 1) 



Z Z a (o-j 1 ,i-i> • 1 

3 X =0 3 2 = 0 

3a + 3 2 ^1 



= a, 



0 1 
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1 1 



c, ’ ,) = Z Z a(i_ ii> i '32 ) 

jj=0 j 2 = 0 
3 1 + : 2 ^ i 



2 ot oi a io 



( 2 , 0 ) = 



1 1 

j 1 = o j 2 =o 

V 3 2 * 2 



(2-j 1 , 0- j 2 ) = 



a 20 + a !0 



'( 2,0 



2 1 

I E 

3i=° 3 2 =0 

3 i + 3 2 i 3 



a (2-j 19 l-j 2 ) b (j 15 j 2 ) 



a iO b (l,l) + a 01 b ( 2 ,0) + a 20 b (0,l) 



2a i0 a 01 + 2a 01 a 20 



These coefficients are, of course, the same as in eq (4.12a). 

The computations for the b^ 2 ^ coefficient can be dis- 
played graphically in the following manner: if the coeffici- 
ents of the 2-variable characteristic equation of L(z 1 ,z 2 ) 
are arranged in matrix form, as follows, 

OtQ 2 * * • 

QL\ z * * * 

Oi-2 2 • • • 
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1 


a oi 


a io 


an 


a 20 


a 2 i 


• 


• 



and if we arrange the b.. . .in matrix form, 

^ 3 i ? 3 2 ' 



1 — 1 


b (o, 1) 


b (°, 2 ) 


( 1 , 0 ) 


b (l , 1 ) 


b (l , 2 ) 


(2, 0) 


b (2 , 1 ) 


b (2,2> 



then the equation defining b^ 2 ^ , 



..... ■ I I 

VO j 2 =o 

3 , + 1 * 3 



a (2-jj,l-j 2 ) b (j 1 ,j 2 ) 



can be visualized as the rotation of the first matrix by 
180 ° and the superimposition on the second matrix, such that 
a 00 = 1 is above b^ 2 1 ^ . By summing the products of stacked 
coefficients, as shown in Figure 4 , the result is: 



b (2,i) ’ a io b (i,i) + a oi b (2,o) + a 20 b (o,i) 



+ a 2 i + ai i b / x 

i i , o ) 



3 . Nonrecursive Formula 



To compute the term t’(j) using the recursive formula 
requires the computation of all previous terms, i.e., b^T^, 
where 3 < I. To circumvent this disadvantage, a nonrecursive 
formula is derived in the following paragraphs in which 
b^T^ is uniquely computed using the set of coefficients in 
Q(z ) . 
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FIG. 4: GRAPHICAL INTERPRETATION OF EQ (4.21) FOR THE 
b (2 1) COMPUTATION OF EXAMPLE 4.2 
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Equation (4.15) can be rewritten as 



L(z) 



00 



I 




Q*(z) 



Then, if we define 

M^(z) = Q^Cz) 
eq (4.23) becomes 



L(z) 



Y, H t Cz) 
1 = 0 



(4. 23) 



(4.24) 



(4.25) 



and taking the (i 15 ...,i^) derivative with respect to 
(Zj,...,Zjj) of both sides in eq (4.25) leads to 



L ( , . . . , l N ) ( z ) _ V’ ( ij , . . . , i N ) ) 

£=0 



(4.26) 



where the notation of eq (4.10) is used. 

It is concluded that the Taylor coefficients for L(z) 
expansion can be computed whenever the term M » ' • • > ( z ) 
can be calculated. 



Theorem 4.2 : 

If R(x) = r, (x)...r (x), where r-(x) is a function of 
1 m 1 

one variable, then as shown in Chrystal [54] 
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U) 

R(x) = 






(x) ... r (x) 

m 



l, + ... 1 = £ 
1 m 



(4.27) 



Theorem 4.2 applied to eq (4.24) results in: 



M, 



(ij,o, . . . ,o)(-) = V 1 / X| V 

Z_ j yl io’ • • • >^!jy 



( 3io’° ’ • ' • ,0 \z) 



3io + - 



. . • • • ’ 0) (z) 



if Theorem 4.2 is applied again to eq (4.28) 



(4.28) 



£ 



> 0 ) ( 5 ) = 



z z 



• I • • • I • * 

W- " * ^1 1 1 3 20 + * • * H 2 2 



(. \ )(. *• . ) 

\ J IQ 5 * * * 5 ^ 1&/ \ ^20’ * * * 5 ^ 2 JJ, ' 



q( jl0>320»° ’ • • ' 0) (z)-Q ( jl£’3 2 £ 5 0, . . • o y z ) 



(4.29) 
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Then by induction, it follows directly: 














(4.30) 

Combining the coefficient term of eq (4.9) and eqs (4.26) 
and (4.30): 








. q(j 10 > • • • > j ^ 0 \ # # Q ( ji2,i • • • j M 

(4.31) 
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Applying (4.16), (4.17) to 




( 0 < ij < °° ) 



a( 0 ) 



0 



the coefficient term a 



(I) 



is computed as 



a 



(i) 




Q 



(ii,. 




Similarly, for eq (4.20), 



(4.32) 



. L(z) 




( O^ij^ 00 ) 



where 



b 



(i) 






(4.33) 



(4.34) 



But by eq (4.31), b^r^ equals 



CO 




( 






1 

^ No ! • • • 3 N£ : 




• Q 



(3 







n ( 3 iji » • • • » 3 n a > 

Q (0) 



(4.35) 



If eq (4-32) is applied to eq (4.35), then the result is 
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summarized as follows: 



Theorem 4.3: 



If the Taylor expansion of Q(z) is 



Q(z) 



W. v> zl ; 



a (0) = 0 



( 0£iv°°) 



and if the Taylor expansion of L = 1/1-Q is 



(4.36) 



L(z) = 




b (i ) z- 



(4.37) 



( O^i^oo) 



then for all i 



b d) 



■I I 



£=0 



o + * ' ‘ 



I 



j No + " * + ^N£“ 1 N 



• a 



( jio»» • • » jjg„) (ji£ v • • j^nj?? 



Using vector index notation, eq (4.38) becomes 
) 



V 







a (j 0 )--- a (j,) 



• , f 



= 1 



(4.38) 



(4.39) 



where the multiple summation is performed over all ^-tuples 
of vector indices (j Q ,...,j^) which sum to the vector index 
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i. The application of Theorem 4.3 is illustrated in the 
following examples: 



Example 4.3 : 

For the two-dimensional transfer function 



L( z , z 



] 

1 " «, 0 z i - a oi z 2 - « 2 o z i" 



(4.40) 



which was also used in the previous examples, the first two 
coefficients of the corresponding Taylor expansion are com- 
puted as follows : 



(i,o) 



oo 

t . I. . I 

^ 0 3i o + • • • + 3l£ ^2 0 + • * • H 



‘32 i 



= 0 



a (3 10 >3 20 > 



where the only nonzero coefficients are 



(i,o) = a io 



(o,i) a oi 



(2,0) = a 20 



(4.40a) 

(4.40b) 

(4.40c) 



For £ = 0 the indices jio and j 2 0 Sum to ii : 1 an< 3 i 2 - 0 

whenever j 10 = 1 and j 20 = 0. 

For £ = 1, the indices sum to ij = 1 and i 2 = 0 whenever 
^3io»3iP and (j 20 » 32 1) belong to the following set of ordered 



pairs : 



Jj 10 .j M ) = {(0,1); (1,0) J- and (j 2# ,j 21 ) = {(0,0)| 

i.e., 

{inf 0 ’ Ojfl}, {j 10 =1 »3n =0 }» and {j 2 =0,j 2 =0)| 

For £ = 2 , following the same line of argument, 

(jio.jii.ji 2> = {(0,0,1) ; (0,1,0) ; (1,0, 0)| and 

Now, if the indices for £ = 0,1,2,..., and the coefficients 
of eq (4.40) are substituted in 



b O,o) 





a 



( j 1 0’ O 2 



q) 



( £ = 0 ) 



X! X! a( 3i0’j 20 ) a (in.j 2 i ) 

OlO + jll O 2 0^" 0 2 1 

= 1 =0 



+ 




jl0 + 0ll + jl2 
= 1 




( jl 0’ O2 <P 



02O + 021 + 3 2 2 
= 0 



a ( 3l 1» j 2 l) 3 ( j 1 2’ O 2 2 ) 



(£ = 2 ) 



+ . . • 



then the equation on the following page is obtained: 
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b (l,o) a (l>o) 

+ a (o, o) a (i, 0) + a (i, o ) a (o, 0 ) + 

+ a (o, 0 ) a (o, o) a (i, 0 ) + a (o,o) a (i,o) a (o,o) + 

+ a ( J, 0 ) a (o, o) a (o,o) 

+ . . • 

= ai 0 , since a ( 0 ^ = 0 . 

Similarly , 

b (o, l ) = a 0 1 



Example 4.4: 



The term, b^ 5 ^ is computed using Theorem 4.3 for the 
transfer function L(z l5 z 2 ) , which is defined in examples 4.3 
as follows : 

For N = 2, i a = 5, and i 2 = 1 , equation (4.38) becomes: 



b (s,i) 





I 



a( 3 10 O2 0 ) ‘ ' * a( 3i ,>j 2 > 



I2 o + * * * 
= 1 



(4.41) 



or 

b (s,‘) = a (s,i) + 



U = 0) 
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I .1 

V^n 5 ^2 o + i l" 1 



a (3 10 > 3 2 0 ) ^n’W (£ = 1) 



I I 



3io + 3,, 
+ 3 12 =5 



a( 3 ln ,3, n ) a (3,,,3, 1 ) a (j,, 5 j,,) 



( £ = 2) 



j 2 o + 3 2 i 

+ 3 22 =1 



•I Z • 

^10 + ^ll + ^12 3 2 0 + 3 2 l + ^2 2 

+ 3 i 3 = 5 + j 23 =1 



( 3 10 4o } a( ^12’3 22 ) • 



( 3i 35 j 2 3 ) 



( £ = 3) 



•£ I 



( 3 10 ’W a( ^11^2 1 ) a( 3l2^22 } 



^ 1 0 + ^ 1 1 + ^ 1 2 ^2 0 + \ 1 + ^2 2 
i +i =5 i -H =1 

J 13 J 14 J 2 3 J 24 ^ 



a (j 13 4 3 ) 



( £ = 4) 



* X! X! a<i >« ,j2 * > a<j n ,j 2> ) a<j .2’i 2 2 ) 



■^ 10 + ^ ll + ^12 ^2 0 + ^2 1 + ^22 

^U^IS 3 2 +3 2 +3 25 
= 5 =1 



a ( 3 1 4 , 3 24 > a (3 15 ,3 2; ) 



( £ = 5) 
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+ ... 

The only nonzero coefficients in the characteristic equation 
of L(Zj,z 2 ) are : 



(l, 0 ) 


P 

O 


(4.42) 


(0,1) = 


a 

01 


(4.43) 


(2, 0) = 


a 

20 


(4.44) 



It follows from eq(4.41) that the sets of indices 
|ji<p an ^ | j^o’* * * » 3 2 £, j 

must add up to 5 and 1, and the set of ordered pairs: 

■| j 1 0’ ^2 ’ * • • ^ 3 1 £_ > j 2 £ ) | 

must correspond to nonzero coefficients, i.e., as listed in 
eqs (4.42) - (4.44) . 

Applying this reasoning to term ( i = 1) the two sets of 
indices are obtained. For 

{(0,5); (1,4); (2,3); (3,2); (4,1); (5,0)} 

and for 

( j 20 » j 2 P » {(0,1);(1,0)| 

It can easily be verified that the products of coefficients, 
which is defined by term ( Z = 1), becomes 
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I I 



< □ 105 320> < jll» D2 l) - 



iio+^f 5 


^2 0 + 3 2 1 1 










a ( °, 0 ) 


a (5,l ) 


+ a / x 

( 1 , 0 ) 


a (4, 1) 


+ 


a (2 , 0 ) 


a (3,1 ) 


+ a , x 

( 3 , 0 ) 


a (2, l) 


+ a (4,0) 


a (l , 1 ) 


+ 


a (5 , 0 ) 


a (o, 1 ) 


+ a (o, 1) 


3 ( 5, 0 ) 


+ 3 ( 1 5 1 ) 


a (4 , 0 ) 


+ 


a ( 2 5 1 ) 


3 (3 , 0 ) 


+ a ( 3 5 1) 


a (2,°) 


+ a (4,l) 


a (i, 0 ) 


+ 


a (5, 1) 


a (o,o ) 



Following the same reasoning, the term ( i = 2) becomes 

Yl X^ a( 3 10 O 20 ) a( 3 11? j 21 ) a (3i?i 22 ) = 0 

-^10 + ^ll + ^12 ^20 + ^21 + ^22 

= 5 =1 

The two sets of indices for the term ( i = 3), which add 
to 5 and 1, as well as corresponding to nonzero coefficients 
in eqs (4.42) - (4.44), lead to 



X! X! a( ^o’3 20 ) a (: 1 pj 21 ) a (: 12 >j 22 > a (3i3^23> = 

-^10 + ^11 + ^12 ^20 + 3 2 l + 3 22 

□l 3 = 5 323 = 1 

t 

= a (2,0) a (2,0 ) a O,0 ) a (0,l ) + a (2,0 ) a (l, 0) a (0, 1) a (2,0) + 



+ a (i,o) a ( o, i ) a ( 2 ,o) a ( 2 ,o) + a (o, i) a (2, o) a (i, o) a ( 2 , o) + 
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+ a (2,0) a (l,0) a (2,0) a (0,l) + a (l,0) a (2,0) a (0 } 1) a (2,0) + 



+ a ( 2, 0 ) a (0 , 1 ) a (2, 0 ) a ( l , 0 ) + a ( 0 , 1 ) a ( l ,0 ) a (2, 0 ) a ( 2,0 ) + 



+ a ( 1, 0 ) a ( 2 , 0) a (2, 0) a (0, l) + a (2,0) a (2,0) a ( 0 ,l) a O,0 ) + 



+ a (2,0) a (0 ,1) a (l,0) a (2,0 ) + a (2,0) a (2,0) a (0,l) a (l,0) 



1 2 ot 20 CL jg Otg i 



Similarly, the term ( £ = 4) is 



.1 




2l 0 + ^l l + ^l 2 


3 2 0 + ^2 1 ^2 2 


i +i =5 

J 13 J 1 4 


^2 3 + ^2 4 - ^ 


— 20 Ct j g 


a 0 1 a 20 


and the term 


( £ = 5 ) becomes 


.1. 


X! a( to^20 > 


3l0 + il 1 + 3l 2 


^20 + %. l + ^22 


ii3 + iur 5 


^2 3 + ^2 



a , . 



a , . 



k a, . 



( to^2 0 )a( 3ll^ 21 ) a (j 12 ,j 22 ) a ( 13 ,j 23 ) 



a (3 1 ^ 24 ) a (j 15 ^2 5 ) 



6 a !0 5 a 0! 



It can be shown that the terms for l > 5 are identically 
zero. 
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Thus : 




These examples show that the application of Theorem 4.3 
is straightforward, but not trivial. The computational 
method allows the determination of some b without requiring 
the knowledge of any previous b's. The following section 
develops a nonrecursive combinatorial formula which has 
several advantages as compared with the above outlined tech- 
nique. 

4 . Nonrecursive Combinatorial Formula 

Following [53] and using the vector notation defined 
in Appendix A, i.e.. 



(j ) = (jj>- • • , j N ) 



it is defined that two Jl-tuples 



( V • • • > j and 




are equivalent, i.e.. 



< i 0 » • • ♦ » j £ § ( j 




whenever any given vector index j appears the same number of 
times in each £-tuple. Therefore, if 



then it can be easily shown that 

(Do ) ^o ) ( 3 £ ) 

Also given an £-tuple of vector indices for each 

vector index j , tT defines the number of times the vector 
index j repeats among It follows directly that 



3. , • N 

OJ 



( V 



n 

j 



lz 

L (j) 



(4.45) 



For example, consider the term 



where 

then 



a (jo) a (j*> = a O,o) a (o, i) a (2,of (i,o) a (i,o) 
j = { O, 0 ) ; (0,1); (2,0)} 

tT = {3 ; 1 ; 1 } and a (Jo) ...a (J } = a ( 1 , 0 ) a ( o , 1 ) a ( 2 } o ) 



It is also observed that the sum of tT is equal to i + 1, 



1 . e . 



I t 5 2 [m] 

j 



(4.46) 



If there are (s or ..,s^) distinct vector indices among the 

i-tuple, (j 0 ,...,j^), then t- i 0 for n = l,...,k. Since 

s n 

another £,-tuple ( , j £) is equivalent to( j 0 , ... ,3 ^) if 



and only if each index (s o ,...,s^) appears in( , . . . , j £) 

the same number of times t- , ... ,t- it appears in (j 0 ,...,j 0 ) 

s o s k x 
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it follows from combinatorial theory that the number of 
H-tuples equivalent to (] 0 ) equals the multinomial co- 
efficient, 




(4.47) 



Note : 

S^ = (ii,...,i^) is defined with k = 0,1,... for each 

distinct index (i , ...,i^) which appears in the coefficients 

of Q(z) . t- is a positive integer. 
s k 

Combining these results, the following theorem is obtained: 
Theorem 4.4 : 

If the Taylor expansion of Q(z) is 



QCzl = >> a (I) 



T . Z 



•’ a (B) = 0 



( 0 < ij <«>) 



and the Taylor expansion of L = 1/1-Q is 




l 



z 



( 0 < i^ <°°) 
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where 





must be satisfied. 



(4.48b) 



Example 4.5 : 

For the two-dimensional transfer function used in the 
previous examples, the first two coefficients of the corre- 
sponding Taylor series expansion are computed as follows: 




(o, i) +t ( 2 , oA t (i , o) t (o, 1 ) 

(O.DV^O )] a(1 ’ # ) a < # ’ 1) 



t (2, 0) 
a (2, 0 ) 



where t (1 ^ 0) (1,0) + t (0 ^ (0,1) + t (2 0) (2,0) = (i a ,i 2 ) 
For ij = 1 and i 2 = 0 



P.,.) (1 ’ 0) + V.) (0 ' 1) + Pa,.) < 2 ’°> = 



whenever t (M) =1 , t (D>1) = 0 , t (2>t) = 0 



Therefore, b 



5 “(i,o) = a (i,o) " ^10 



a. 



Similarly, b (0 jj 



a 



01 
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Example 4.6 : 

For the transfer function , L(z), as in example 4.3, 
compute b^ 5 ^ . The only nonzero coefficients in Q(z) are 



Thus , 



a d,o) - 
a (°, 1 ) = 

a (2,0) 

S 0 = (1,0) 
ij = (0,1) 



s 2 = (2,0) 



a 



a 



10 



o i 



a 



20 



and t (lj0) (1,0) + t (2j 
eq C4.48b), whenever 


o) <2,°) 


+ t C,i 


) (0,1) 


= (5, 


1) from 


t O,o) = 


1 


t (2,0) 


2 


t (o,i) 


= 1, 


or 


^ (l , 0 ) 


3 


^( 2 , 0 ) 


1 


t (o,i ) 


= 1, 


or 


t (l,0 ) = 


5 


t (2,0) 


0 


t (0,D 


= 1. 





Therefore, substituting into (4.48a), 



( 5 , 




Q 10 a 01 a 20 



3 + 1 + 1 



3,1,1 



a i0 a 0 1 a 2 0 + 




a io “oi 



12 a io a oi a 20 2 + 20 a i0 3ot tti a 20 + 6 a 10 a 



01 
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5 . Comparison of results 



The four computational schemes presented in Part C 

are summarized as follows: 

The derivative method is very difficult to apply 

since the computation of the coefficient of the Taylor 

~th 

expansion of L(z) requires taking the (ij,...i^) partial 
derivatives of L(z) and evaluate the result at zero. This 
technique cannot be easily implemented on a digital computer. 

The recursive method requires for the computation of 
b(I) knowledge of all b^t^, where j <i. This technique 

is best suited for computer implementation, especially when 
computing a finite set of coefficients, i.e., {b^:jj} > where 
0 < i^ < L. It will be used extensively in Chapter V. 

The nonrecursive and nonrecursive combinatorial 
formulas relate the coefficient b^j^ explicitly to the co- 
efficients of the N-dimensional characteristic equation. 

To compute b ,, . x for the transfer function, 

^ K 5 3 > 



L(z, , z „ ) = 



1 ’ 2 • 



1 - a io z i - 



a 01 Z 2 a 20 z l 



Using the nonrecursive technique of Theorem 4.3 requires 




m= 0 



(k-m+j ) ! 

JT ml (k-2m) ! 
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multiplications in addition to identifying the same number 
of index sets, where [X] means the largest integer < X. 

When applying the nonrecursive combinatorial method 
of Theorem 4.4, only + 1 multiplications must be per- 
formed in addition to identifying number of t's and 

multinomial coefficients, which are tabulated in [55]. The 
following table further elucidates the comparison of the 
number of multiplications to obtain b^ ^ for the above 
L(z ,z ) : 



<k,j> th 
k = 


term 

j = 


Number of computations 
Th 4. 3 


Th 4.4 


0 


0 


1 


1 


1 


1 


2 


1 


2 


2 


9 


2 


3 


3 


40 


2 


4 


4 


190 


3 


5 


5 


924 


3 


10 


10 


3045185 


6 


15 


15 


1.1473 E + 10 


8 


20 


20 


4.5734 E + 13 


11 


25 


25 


1.881 E + 17 


13 


30 


30 


7.89227 E + 20 


16 



The tremendous savings in computations when apply- 
ing the nonrecursive combinatorial method of Theorem 4.4 are 
obvious . 

There is one important application of Theorem 4.4 
which will be used in Chapter VII to establish a mathematical 
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formulation for the propagation of transfer function coef- 
ficients in unit sample responses in any dimension and any 
order. These rules are then. used to develop in Chapter VIII 
a transfer function extraction algorithm, which in turn is 
used in Chapter IX to develop time domain design 
techniques. Therefore , the high speed characteristic 
which evolves from the few multiplications is exceptional 
but the mathematical form is of greater importance for the 
subsequent chapters of this thesis. 

D. UNIT SAMPLE RESPONSE 

By comparison of eqs (4.13) and (4.19), the relationship 
H(z) = L(z -1 ) is easily identified. The following theorem 
establishes the relationship of the unit sample response 
entry h(i) and b , where h(I) is also the I th coefficient 
of the Taylor expression of H(z). 



Theorem 4.5: 



If the Taylor expansion of Q(z) if 



Q( z ) = 




a (i) Z ’ 



; a (0) = 0 



, 1 
^ 0 <i_^ <<»^ 



(4.49) 



and if the Taylor expansion of L(z) = - 



1 - Q(z) 



is 



L(5) = 



b d> zl 



(4.50) 



/ X 

(o<i.<~) 
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1 



in the 



then for the Taylor expansion of H(z) 
form 



H(z) 






-l 

z 



1 - Q(z 




(4.51) 



each and every b^T^ = h(i) for all i, and the coefficients 
h(i) are related to a^^ as stated in eq (4.17) and in 
Theorems 4.1, 4.3, and 4.4. 



Note that the notation h(i), instead of h^^j, is the 
commonly accepted one to designate the i^ coefficient of 
the unit sample response. 

Table 4 summarizes the h(i) versus a^^ relationships. 

Equation (4.51) can be written for the two-dimensional 
case as matrix product: 



where 



= (1 


z ' 1 z ~ 2 
Zj Zj 


. ) A (1 


z 2 -1 z 2 2 . 


h ( 0 , 0 ) 


h(0 ,1) 


h ( 0 , 2 ) 


h ( 0 , 3 ) 


h( 1 , 0 ) 


h(l,l) 


h ( 1 , 2 ) 


h(l, 3) 


h( 2 , 0 ) 


h( 2 , 1 ) 


h( 2 , 2 ) 


h ( 2 , 3 ) 


h ( 3 , 0 ) 


h ( 3 , 1 ) 


h ( 3 , 2 ) 


h ( 3 , 3 ) 



and h ( 0 , 0 ) = 1. 



) 



t 
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TABLE 4: FORMULAS FOR COMPUTING UNIT SAMPLE RESPONSE 

ENTRIES h(n) FROM A GIVEN TRANSFER FUNCTION H(z) 
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v - app lication of series method for the derivation 
OF STABILITY C ONDITIONS IN TWO- DIMENSIONS 



A. INTRODUCTION 

The formulation developed in Chapter IV is applied to 
compute the unit sample response matrix {h(n ls n 2 )} of the 
third and fourth order uncoupled, as well as second order 
coupled transfer functions. Conditions for absolute 
convergence of matrix {h(n 1 ,n 2 )}, i.e., 

CO CO 




which are proven in Appendix C to be identical to BIBO 
stability conditions, are found in two steps: 

(a) for the partial infinite sum of absolute entries 
in each column (row), and 

(b) for the infinite sum of column (row) sums, where 
the sequence of summation steps corresponds to Cauchy's 
method one (two) respectively. It is shown in Appendix K, 
that the summation of successive diagonals, which is Cauchy's 
method three, leads to the same result. 

For the third and fourth order uncoupled transfer function 
the latter method is used in conjunction with binomial theory, 
resulting in a sequence of geometric series for which 
convergency conditions are known. In order to apply the 
binomial theorem m forming the sum of successive diagonals 
{ |h(ni ,n 2 ) | } , the absolute signs must be removed. This is 
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done, except for the trivial case, by identifying each entry 
h(n ls n 2 ) to represent an Ultraspherical polynomial. By 
using its known characteristics, necessary and sufficient, 
as well as sufficient stability conditions, are derived 
(see tables 5.1a, 5.1b). 

It is also shown that in the second order coupled case, 
Cauchy's methods one or two in conjunction with derivative 
operator theory results in a sequence of geometric series 
for which necessary conditions for column convergence and 
necessary and sufficient conditions for column sum convergence 
are identified (see table 5.1c). 

B. UNCOUPLED FILTER 

Definition: Uncoupled Filter 

A two-dimensional digital transfer function is uncoupled, 
if its two-dimensional characteristic equation has no cross 



n 2 are both non zero. 

The general uncoupled two-dimensional transfer function 
can be written as 



terms, i.e., terms of the form a 



nin 2 




where n i and 



H( z i ,z 2 ) 



1 



(5.2) 



g( z i ) - h( z 2 ) 



00 




(5.3) 



p = 0 
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Transfer Function H(z,,z,) Case Stability Conditions 
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1 . Third Order Uncoupled Transfer Function 

-1 -2 

For g(zi) = 1 - aiozi - a 2 oZi 



(5.4) 



h(Z2) = CX01Z2 1 



(5.5) 



eqs (5.2) and (5.3) become 

H( z i , z 2 ) = ^2 IT (5.6) 

1 - aiozi - a 2 oZi - a 0 iZ2 

00 

= (1 - aioZi^ - a 2 oZi^) auz 2 ^ 

( 5 . 7 ) 

p=0 

Equation (5.6) is realized in direct form as shown in Fig 
5.1. Note that the direct forms are identical whenever the 
numerator polynomial is equal to unity. Equation (5.7) can 
be used to implement the third order uncoupled recursive 
filter in hybrid form, i.e., non-recursive structure with 
recursively realized components as shown in Fig 5.2. 
a. Unit Sample Response 

The planar unit sample response A of the two- 
dimensional third order transfer function (5.6) is listed 
in Table 5.2. The details of computations, which are 
performed using the recursive formula of Theorem 4.1 are 
contained in Appendix I, Section A. 



106 




a 

01 



Fig. 5.1= DIRECT FORM REALIZATION OF THIRD ORDER 
UNCOUPLED FILTER IN TWO-DIMENSIONS. 
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I 

I 





yCn,,^) 
— D> 



Fig. 5.2= HYBRID STRUCTURE OF THIRD - 
ORDER UNCOUPLED FILTER IN 
TWO-DIMENSIONS. 
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b. Summation of Absolute_Entries of 
Unit Sample Response A 

The sum of all absolute entries of column 0 

becomes : 



S o 




|h(n x ,0) | 



V 

m=0 



, \ . ni-2m m 

(ni-m) ! ajo a 2 o 



m 



! (n i-2m) ! 



(5.8) 



Similarly, 



S i 



I 

n j = 0 



[n i /2 ] 



m=0 



(n t - m+1) ! a7o~ 2m an i 
m! (n !-2m) ! 



m 

£Lz jl 



(5.9) 



and for the 



th 

j 



column 



OO 


[n,/2] 


« V 


t — , , . N , n.-2m j m 

\ ( n,-m +] ) l aio an i a,n 


L 


/ j ! m! (n,-2m) ! 


n i = 0 


m=0 



(5.10) 



The sum of all absolute entries of the unit pulse response 



oo oo 




n 2 = 0 n x = 0 



is now easily identified as: 



oo [n x / 2 ] 



= V V 


\ / , \ . n i “ 2 in n m 

\ C n j +n 9 -m) . aiQ ao i_ ot?_n_ 


z_ 2_ 


/ n 2 ! m! (nj-2m) ! 



n 2 -0 n x = 0 m=0 



(5.11) 
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c. Necessary and Sufficient Stability Conditions 

2 

for CX 20 ^ - — jp- using Ultraspherical 

Polynomial Method 

It is proven in this section that the two- 
dimensional third order transfer function (5.6) is BIBO 
stable if and only if 



0t2 0 < 1 - | ct 1 0 | - | cx 0 1 | 



(5.12) 



2 

whenever a .20 ^ by inspection of eq (5.11) first for 

the case (X 20 0 and secondly by utilizing the known 

characteristics of Ultraspherical polynomials. 

(1) Stability Conditions for 0620 >. 0 . It is 
proven in Appendix C that a two-dimensional system is BIBO 
stable if and only if 

OO 00 

|h(m > n 2 ) | < » (5.13) 

n 2=0 ni=0 




00 00 



Equation (5.11) can be rewritten as: 
[n,/2] 



S = 



) ) ni n 2 ) (m+n 2 -m) ! / a 20 \ 

L - — Z_ 10 01 n 2 ! m! (ni-2m) ! ( 2 ) 

n 2 =0 m=0 1 m=0 X ai ° ' 



m 



(5.14) 



111 



which for a 



> 0, becomes 



2 o 



oo oo 



s = 




[ni/?] 




(rti+n2-m) I 
ml (n!-2m)! 




n 2 -0 ni-0 



m=0 



(5.15a) 



since an and ao i are real numbers. 

The summation of all entries of <A> , 
where < > is defined as taking the absolute value of each 
entry, i.e., the (n l5 n 2 ) L term; using eq (5.15a) is 



|h(ni ,n 2 ) | 




(ni+n 2 -m)! 
n 2 ! m! (n 1 -2m) ! 



m=0 




(5.15b) 



is performed following Cauchy's third method of summing all 
entries of an infinite by infinite matrix. (See Appendix K) . 
With the detailed computation contained in Appendix I , it is 
noted that the sum of successive diagonals of < A > is 
greatly simplified by application of the binomial theorem 
and derivative operator methods . 

The summation of all entries of A 
converges or equivalently the third order filter transfer 
function H(zj,z 2 ) is stable if and only if 

CI 20 < 1 — | ot i o | - | oc o i | for ot 2 o 0 • (5.16) 
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Notice that the validity of eq (5.16) is easily checked by 
inspecting the unit sample response matrix A as shown in 
Table 5.2. 

2 

(2) Stability Conditions for 0 > ct 2 o > “ 

It is shown in this paragraph that eq (5.14) can be rewritten 
for a 2 o < 0 such that the summation inside the absolute signs 
can be identified as an Ultraspherical polynomial. In order 
to simplify the summation of all absolute entries of A by 
application of binomial and geometric series theory it is 
necessary to remove the absolute signs, which can be done 
utilizing the known characteristics of the Ultraspherical 
polynomials . 

For a 2 o < 0 , eq (5.11) can be rewritten (5*17) 
in the following form: 



II 



n 2 =0 ni =0 



tm/2] 



|(/-a 20 ) ni oft? 



I 



m=U 



(- 1 ) 



m (ni+n 2 -m) ! 
m! n 2 ! (ni- 2 m)l 



ai o 

/“ 0 t 2 ( 



nj- 2 m 



(5.18) 



But the summation inside the absolute sign is identified 
using the theory presented in Appendix D, as Ultraspherical 
polynomial , i . e . , 



(n 2 +l) 

C 

Hi 



Olio 

2 /-a z o 



[m /23 





(ni+n 2 -m) ! 
”n 2 ! (ni- 2 m) 



i 




m=0 



ni-2m 

(5.19) 
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It is proven in Appendix D that the Ultraspherical polynomial 



(n 2 +l) 



ni 



/ ai 

\2/Ta 



2 0 



is positive, 



whenever 



ai ° > 1 , 
2 /-a 2 o / 



(5.20) 



which implies aio >0 . Therefore, eq (5.18) can be 



rewritten as 

OO 00 

..II 

n 2 =0 ni=0 





(n 2 +l) 

C 

n l 



( j (5.21) 

\ V-a 20 / 



Also, for 



( . -gj£ ) <-l 

\ 2 /- a 2 o / 



it is proven in Appendix D that 



(n 2 +.l) 





(5.22) 



is positive/negative whenever n t is even/odd, respectively. 
Thus in this case, the absolute value of the Ultraspherical 
polynomial becomes equal to 



(5.23) 



(n 2 + l) 

C 

n i 



O I 0 



,2/T 



0.2 0 



(-I) 111 



(n 2 +l) 

C 

n i 



a i o 



„ 2 /-a 2 o 



It is also proven in Appendix D that the right hand side of 

*. 

eq (5.23) is equal to 



(n 2 +l) 




-c* i o 
2 /-a 2 o 



(5.24) 



Condition 



o 1 o 



< -1 



1 2 /-a 2 o 



(5.25) 



implies ai 0 <0 ♦ Thus 



(n 2 +l) / \ (n 2 +l) 

; *JJ1. . = C -1SLU. 

ni \2/-a 2 o/ Rl V 2 /- a 2 o 



and 



00 00 



■ II 

n 2 =0 nj=0 



-n, , ,n 2 „ (n2+1) / |q, q 1 



(^5To) 1 la,,!" 2 c 



ni 



1-J- 



(5.26) 



a 2 o 
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Conditions (5.17), (5.20) and (5.25) are now combined as 



0 > CX2 o > 



a i o 
4 



(5.27) 



and similarly, eqs (5.21) and (5.26) to 



S = 



00 00 

II 



( /-ct 2 o) ni 




(n 2 +l) 

r> 

ni 



n 2 =0 n!=0 




which by eq (5.19) is equal to 



00 [n i / 2 ] 




n 2 =0 n x =0 m=0 



(ni+n 2 -m) ! 
n 2 l m! (ni-2m) ! 




(5.28) 



which is identical to eq (5.15). Therefore, it is concluded 
that H(z i , z 2 ) as defined in eq (5.6) is stable if and only if 

a 2 o < 1 - | ot i o | - | ot o l | (5.29) 



for 



0 > a.20 > 



aio 

4 



(5.30) 
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d. Sufficient Stability Conditions for 0120 < - 
It is proven in this section that the third 
order uncoupled transfer function is stable if 



aio 

4 



| oto 1 | < (1 - /-a 2 0 ) 2 



for 



a 2 0 < 




( 5 . 31 ) 



by using the known characteristics of Ultraspherical 
polynomial. For 



a 2 0 




5 



which is the condition for complex roots of the one- 
dimensional second order transfer function: 

H( z i ) = — 

1 - aizT 1 - a 2 z j 2 

it is shown in Appendix D that the Ultraspherical polynomial 



(n2+l) 





t 



has ni real and distinct roots. This implies that the sum 
of all absolute entries of "K must be performed without 
removal of the absolute signs enclosing the Ultraspherical 
polynomial part in: 
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00 



00 




n 2 =0 n^O 



n 2 



(n 2 +l) 




a io 

2 /-a 2 o 




(5.32) 

/ 



If we apply Cauchy's first method of summing all entries of 
an infinite by infinite matrix (See Appendix K) then the 
sum of all entries in column 0 becomes 



00 



So 




n x 

( / - 0t2 0 ^ 





(5.33) 



It is obvious from Fig D.l of Appendix D that there exists 
'-—Oh 11 . \ -1 -n o(l) 



where all C'" w are of same sign. 

n i & 



no range of 

\ 2/— a 2 o 

Therefore, an approximation from Ref [55] is used, i.e., 



Thus 



( n*+l) 



n. 



-QLul 



^2 /—ct 2 o 



n 1 +2n 2 +l 



n, 



(5.34) 



s t < Y. 



n ! = 0 



( / Ct 2 0 ^ 



n i / n j + 1 



n, 



/= 



d y~a 



2 0 



.fl 2-0. 



1 - /=5 



(5.35) 



2 0 
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if and only if | /~a 2 o I < 1 




1 

/^aTo 



2 



By similar computations, 



S- 

3 



becomes : 



S. 



< 



(1 




Ct 2 0 ^ 



j 

2 jT 2 



and 

00 



OO 



s = 




s. 

3 




( 5 . 36 ) 



which converges , if and only if 



2 

|aoi| < ( 1 - /-OL20 ) 



q.. e . d . 



e. Improved Sufficient Stability Conditions 
for a 20 < ~ nj -11 - 

It is proven in Appendix I that a two-dimensional 
third order (uncoupled) transfer function is stable if 



TT 

g 

sin 0 • y / ~a 2 0 1 6 2 + l 112 « 2 0 ] [1 ~ > /C a 2 o" 1 

|«oi| TT 

C • e [1 + /=5£7 e ] 



( 5 . 37 ) 
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where 



0 = arc cos 



_2Lul 



2/T 



c = 1.11 



(5.38) 



a 2 o 



and 



<*2 o 



< _ a l o. 

4 



(5.39) 



Equation (5.37) is an exceptionally good approximation for 
the Ultraspherical polynomial in the domain: 



OLLfl- 

2*^577 



< l 



(5.40) 



The resulting condition (5.37) has exhibited, in a large 
number of tests, excellent agreement with experimental 
results . 

The stability conditions for the third order 
transfer function in two-dimensions are summarized as follows 



H(zj,z 2 ) = - 

1 — 01 1 0 z 1 - <*0 1 Z 2 1 “ <*2 0 Z 1 2 

a 2 

represents a stable algorithm, when a 2 o > ~ jJ -0 - if and only 

if 

<*2 0 < i “ l a lo| - | CX 0 1 | 

2 

when ct 2 o < “ if i s sufficient that 

^ 2 

I CX 0 1 I < - / - CX 2 0 ^ 

or alternatively 

sin 0 / r 57o’ [e 2 + ln 2 v c a7o'] [i - v'^aTo’ 8 ] 

l«0 1 | £ TT 

c 0 (1 + a 2 o ) 

where e = arc cos ( — 2 Li_q — j and c = 1.11 

' 2 /“oTo* } 
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2 . 



Fourth-Order Transfer Function 



For g(zj) as defined in eq (5.4) and 



h ( Z 2 ) - Ot q 1 z 2 * t Cl o 2 z 2 2 



equations (5.2), (5.3) become 



H(zj ,z 2 ) 



1 

1 — aioz" 1 “ l Z 2 1 — a 20 zi 2 ~ & 02 z^ 2 



and 



(5.41) 



OO 

[aoizi 1 + a 0 2 z I 2 ] P [1 - aiozT 1 

p=0 





(5.42) 



Equation (5.41) is realized recursively in direct forms 3 
and 4, as shown in Fig 5.3. A hybrid structure realizing 
eq (5.42) is shown in Fig 5.4. 

a. Unit Sample Response 

The unit sample response is computed in Appendix 
I for the fourth-order filter in two dimensions using the 
same approach as in paragraph 2a. The result is listed in 
Tables 5.3 and 5.4. 

b. Summation of Absolute_Entries of 
Unit Sample Response A 

The sum of all absolute entries of Tables 5.3 
and 5.4 will be obtained using Cauchy’s first and second 

1 

method , i . e . , 
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Fig. 5.3 : DIRECT FORM REALIZATION OF FOURTH- 
ORDER (UNCOUPLED) FILTER IN TWO- 
DIMENSIONS. 
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Fig. 5.4 : HYBRID STRUCTURE OF FOURTH ORDER UNCOUPLED FILTER 
IN TWO-DIMENSIONS. 
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UNIT PULSE RESPONSE 2T OP FOURTH ORDER UNCOUPLED FILTER 



n 2 = 


h(k,n 2 ) 


0 


1 — » 
K j. X 

M is 

1 1 


/, . . k-2i i 

(k-i) ! a 2 o ot 2 o 

1! (k-2i ) ! 




i=0 




1 


[k/2] 

z 


( k— i+1) ! ct i o * oc ^ o ot o i 

i! (k-2i ) ! XT' 




1 = 0 




O 


[k/2] 

Z 


(k-i+2) ! a^"o 2i a 2 'o a 0 i , (k— i+1) ! a 2 ~ 0 2i a 0 2 


c. 


1! (k-2i) ! 2 ! i! (k-21)! 1! 


• 


1=0 


• 


• 

• 

j 


l\) 

t 1 


• 

• 

k-21 1 .v. j-2ra m 

aio a 20 (j-mfk-i)! a o i a 02 

i! (k— 21) ! L m! ( j-2m) ! 




i=0 


nr=0 



TABLE 5.4: m = k ENTRIES OP UNIT SAMPLE RESPONSE 

CORRESPONDING TO FOURTH ORDER UNCOUPLED 
TRANSFER FUNCTION 
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I 



Summation by columns: 



00 00 I [n 1 / 2 ] 



S C = 



\ \ V CXlO CX20 T 

L L I L i! w,; L 



m-2i i [n 2 /2] n 2 -2m m 

(ni+n 2 -m-i) !a 0 1 a 0 2 



m! (n 2 -2m) ! 



n 2 =0 ni=0 i=0 



m=0 



( 5 . 44 ) 



Summation by rows, respectively: 



00 00 



|[n 2 /2] 



S = 



I III 



nx=0 n 2 =0 m=0 



n 2 -2m m [ni/2] ni-2i i 

aoi a 02 t — (ni+n 2 -m-i) ! aio a 20 

m!(n 2 -2m)! \ i! (m-2i) ! 



i=0 



( 5 . 45 ) 



(1) Necessary and Sufficient Stability Conditions 

It is proven in Appendix I that the fourth order 
uncoupled transfer function is stable for a 2 o 0 and OC 02 0 , 

i . e . , 



S c = 



II 



ni 



n2 



ctio 



cto 1 



n 2 =0 ni=0 



lm/2\ 

l 



i=0 



Ct2 0 
aio 



[n 2 /2] 



m 



i !(ni-2i) ! 



(ni+n 2 -m-i) ! 
m!(n 2 - 2 m) ! 



ap2 

^gpi 



m=0 



( 5 . 46 ) 



is finite, if and only if 



g 2 p + gp 2 < 1 - | ex 1 0 ! - | g 0 1 1 



( 5 . 47 ) 
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2 

In the case where 0 > ao 2 > ~ and a 20 1 0 the summation 

of columns S c is written in the following form 



[n : /2] 



■III 

n 2 =0 m=0 i=0 



m~2i i 

Ohfi a2 0 _ ( n i-i H 

i! (ni-2i)! V 




m=0 



/• _-> \ m (ni+n 2 -m-i) ! 
m! (n 2 -2m) ! 




But the term in the square bracket defines the Ultra- 

S "t 

spherical polynomial of the (m-i+1) order, i.e., 






[n 2 /2] 

T I 



(- 1 ) 



m (ni+n 2 -m-i) ! / Ooi 



m=0 



m! (n 2 -2m) ! 






ao2 



n 2 -2m (ni-i+1) 

s 



ao i 



(5.48) 



and 



S c 



00 00 

II 

n 2 =0 m=0 



tm/2] 




i=0 



ni-2i i 

cti o a 2 o (ni-i) ! 

i! (ni-2i) ! 




n 2 (m-i+1) 





It is proven in Appendix D that the Ultraspherical 
polynomial as defined by eq (5.48) is positive for 



( — a ° - l - \ >1 (5.49) 

\2 /-ao 2 / 
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and positive/negative for 



Otn i 



2 /*■ a 



< -1 



0 2 



for even/odd n 2 respectively. 



(5.50) 



Following the same reasoning as in the previous section, it 
can be shown that for a 2 o >0 , 



(m-i+1) 

C 

n 2 



( 



2 /^ 77 / 



> o 



whenever 



0 > a 0 2 > 



«ol 

4 



(5.51) 



(5.52) 



Similarly it can be shown that the summation of all entries 
of <A> by rows, i.e.. 



00 00 
ni=0 n 2 =0 



[n 2 /2] 

I 

m=0 



n 2 -2m m 

ap l ap 2 

m! (n 2 -2m) ! 



[ni/2] 




ni-2i i 
(ni+n 2 -m-i)! ai 0 a 20 
il (ni~2i)! 



equals 



for a 2 p < 0 , 





M2] 

1 



n 2 -2m m ni 

ap i ap 2 (n 2 -m) ! ( Aa 2 0 ) 

m! (n 2 -2m) ! 



(n 2 -m+l) 



n 



i 




ni=0 n 2 =0 hfO 



(5.53) 
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and that 



(n 2 -m+l) 

C 

n i 



( _ LslulL \ > o 

V 2 /-a 2 o ' 



for 

0 > a 20 1 . 

In summary, for the conditions 

(5.54) 

(5.55) 



a 2 o > 0 and 0 



ao 2 



> _ 



a 0 i 



ao 2 > 0 and 0 



Ct2 o 



> _ 



Olio 



which are shown in Fig 5.5, the column series S as defined 



by eq (5.44) can be written in the form 



00 00 


foi/2] 


V V m 

/ / l«io| 1 CXo 1 1 


I 


n 2=0 ni=0 


i=0 


[m/2] 


m 


^ — (m+n 2 -m-i)! 


/OlQ 2 \ 


. \ m! (n 2 -2m) ! 


\ Ot 0 1 / 


m=0 






i! (ni-2i)! 



\ 



(5.56) 



which is identical to eq (5.46) and for which the conditions 
for convergence are as stated in eq (5.47). 
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(2) Sufficient Stability Conditions 



For conditions : 



0 > (X2 o i. 



2 

i o 
4 



2 

and 0 > ao 2 ^ - — q~ 



(5.57) 



and for 



^ 2 
a 2 o < - ^ 2 - and/or a 0 2 < - 



(5.58) 



which are also shown in Fig 5.5, the triangle inequality 
is applied to 




n2=0 m=0 i=0 



m-2i i 
aio 0-20 
il (ni-2i)I 



[n 2 /2] 

I 



m=0 



n 2 - 2 m m 
(m+m-m-i)! api ap 2 
m! (n 2 - 2 m) ! 



such that 



S u = 



III 



[n iZ_ 2] lotio | ni (*~T ) [ ^ 23 (ni+n 2 -i^i)!la 01 | ni (-~^) 



i l (n -2i) ! 



m! (n?-2m) ! 



n2=0 ni=0 i=0 



m=0 



(5.59) 



It is easily observed from the derivation leading to 
necessary and sufficient stability conditions, that the 
conditions for convergence of eq (5.59) must be 

| ot 2 0 | + I ot 0 2 | < 1 - | ot 1 p | - | oto 1 | • (5.60) 
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Eq (5.54) 



Eq (5.55) 



Eq (5.57) 



Eq (5.58) 



Fig. 5. 5 1 




GRAPHICAL INTERPRETATION OF 
INEQUALITIES (5.54, 55, 57, 
58). 



131 



C. COUPLED FILTER 

1 . Second Order Coupled (Bilinear) Filter 

The transfer function of the second-order coupled 
filter, which is also referred to as bilinear case in 
Refs [13] and [24] has the form: 



H( zi ,z 2 ) 



1 

1 - cxioZ^* — aoiZ 2 * — a! izl ^2 1 



(5.61) 



Equation (5.61) can be realized in direct form as shown in 

Fig 5.6. It is noted that for a unity numerator the 

realization in Fig 5.6 represents all direct, as well as 

the canonic forms. 

a. Unit Sample Response A 

The unit sample response A , which is defined 

in Section D of the previous chapter, is computed using 

Theorem 4.1. The resulting matrix is listed in Table 5.5, 
't'h 

in which the k entry of column 1, i.e., 

(k+1) an ao 1 + k an ^ aii 



is rewritten as 

1 

1-m m 

a 1 0 aii 

(1-m) ! in! 

m=0 

d k 

where — ( ) 

d an 




d 

d ai 0 



( a 



k+l-m 
1 0 



stands 



for k 



th 



derivative of ( ) 
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Using the same notation, the j”*"* 1 entry of the row can 



be written as 
k 



I 



m=0 



k-m m 
a 1 o aii 

(k-m) ! m! 



a , „ j+k-m^ 

( a ol ) 

d ot o 1 



(5.62) 



b. Sum of Absolute Entries of A 

The sum of absolute entries of A is equal to 



00 00 



s = 



II I 

n 2 =0 ni=0 m=0 



n 2 



n 2 -m m 
ctio otn 
"(n 2 -m) ! m! 



i n 2 



d a?? 



ni+n 2 -m 
«o 1 



(5.63) 



To identify convergency conditions of eq (5.63), both the 
conditions for convergence are found for each particular n 2 
and 0 _< ni <_ 00 , and the sum S to which it converges, as 
well as the conditions for convergence of 

n 2 = 0 




This approach follows Cauchy's second method of matrix 
summation. In particular, the sum of all entries in column 0 
of {h(ni,n 2 )} is shown in the corresponding section C of 
Appendix I to result for aio 0 in 

OO 




ni = 0 
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which converges if ( a 1 o | < 1 to 



So 



1 



1 - a i o 



Si, the sum of entries in column 1, is proven in Appendix K 
to equal, for |a 10 | < 1 , to 



Si 



> 



1 



1 

oti o 



a l o ctp i + aii 
1 - a i o 



cto 1 



Following the same procedure, 



the S j 



becomes 



S. 

1 



> 



1 

1 - Otl 0 



( a i o cxq i + aii + 
V 1 — a i o 




( 5 . 64 ) 



But the sum of all absolute entries of A is defined by 




n 2 = 0 



which is now identified for ) a. i o | < 1 to equal 



I 

a i o L — 

n 2 = 0 



11 2 

/ a i o ap i + aii + „ \ 

\ 1 - a i o 1 1 ao1 / 



( 5 . 65 ) 
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But the right hand side of eq (5.65) converges if and only if 



a i o a o i * a n + 
1 - a i o 



«o 1 



< 1 



or 



dp i + di i 
1 - ai o 



(5.66) 



2 . Discussion 

It will be shown in the next chapter that another 
necessary stability condition for the bilinear case is as 
follows : 



dpi ~ d; i 

1 + ai o 



< 1 



(5.67) 



Moreover it is known [13] that the conditions in eqs (5.66) 
and (5.67) are also sufficient. 
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Fig. 5.6= TWO DIMENSIONAL SECOND ORDER 
COUPLED FILTER IN DIRECT FORM 
REALIZATION. 



TABLE 5.5 ; UNIT PULSE RESPONSE A FOR SECOND-ORDER COUPLED FILTER 



- e» 




I 
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VI. NECESSARY OR SUFFICIENT STABILITY CONDITIONS 



IN N-DIMENSIONS 



A. INTRODUCTION 

It was shown in Chapter V how conditions for absolute 
convergence of all unit sample response entries , which were 
shown in Appendix C to be identical to BIBO stability 
conditions , are derived for the third and fourth order 
uncoupled and second order coupled case. The removal of 
absolute signs enclosing h(n 1 ,n 2 ) in the Taylor expansion 
of a two-dimensional transfer function is possible in general 
by application of the triangle inequality or whenever all 
coefficients in the denominator polynomial are positive. 

This will be shown in the following paragraphs. 

1 . Necessary and Sufficient Stability 

Conditions for Positive Characteristic 
Equation Coefficients in Two-Dimensions 

The summation of <A> entries is easily made for 

the case where all coefficients of the two-dimensional 

characteristic equation are positive. In fact, the sum of 

all entries is obtained by observing that 



H(z i ,z 2 ) 



1 



( 6 . 1 ) 



1 




a 



m i in 2 




(mi+m 2 i 0) 



138 



H(z x ,z 2 ) 



OO 



( 6 . 2 ) 





(mj+m^O) 



The Taylor expansion of H(z 1 ,z 2 ) is defined in eq (4.51) as 



00 oo 

I I 



nj=0 n 2 =0 



h(n x ,n 2 ) z 1 ni z z nz 



(6.3) 



Therefore, each and every h(n 1} n 2 ) can be identified from 

eq (6.2) by computing the binomial expansion and combining 

terms of equal powers in z. and z 2 . Moreover if all a > 0 

m 1 ni2 — 

then by theorem 4.4 h(n ls n 2 ) > 0 and it follows directly 
that 



1 1 

nj = 0 n 2 = 0 



|h(n x ,n 2 ) | 



00 00 

h(n x ,n 2 ) 

n x =0 n 2 = 0 




(6.4) 



It is proven in Appendix C that H(z 1 ,z 2 ) is BIBO stable if 
and only if 



00 oo 

|h(nj ,n 2 ) | < eo (6.5) 

nj=0 n 2 =0 
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or by application of eq (6.4) if and only if 



OO CO 

h(ni ,n 2 ) < 

n i = 0 n 2 = 0 




( 6 . 6 ) 



But the sum of all unit sample response entries is also 

identified from the right hand side of eq (6.3) for 

Zi = z 2 = 1 » which by equating eqs (6.2) and (6.3) equals 

CO CO 

II 

n =0 n =0 



r Mi M 2 



h(n i ,n 2 ) 



£=0 



m^O mj = 0 



°m jm 2 



(6.7) 



Therefore it is stated that H(z 15 z 2 ) 
“m.m, > 0 ’ if “d only if 



is BIBO stable for all 



Mj M 2 



V 

£=0 




(mj+m^O ) 






£ 



< CO 



( 6 . 8 ) 



But necessary and sufficient conditions for convergence of 
eq (6.8) are easily identified as 

Hi M 2 

II ... . > 

■ *» > — in j m ^ 

m 1 =0 m 2 =0 

(m 1 +m 2 ^0) 



1*10 



2 . Sufficient Stability Conditions in Two-Dimensions 



The sum of all <A> entries is by the triangle 
inequality less than or equal to the sum of all entries of 
the unit sample response matrix in which each coefficient 
is replaced by its absolute value, i.e.. 



Mi 



M 2 



I I 






| h(n ! ,n 2 ) | < l_ 

nj = 0 n 2 = 0 £ = 0|_ mi = 0 m 2 = 0 

(m 1 +m 2 ? ! 0 ) 



a. 



m, m 



( 6 . 10 ) 



1 JU 2 



But the right hand side of eq (6.10) converges if and only 
if 



Mi M 2 

/ / la I < 1 , 

L L 1 miiiij 1 ’ 

mi=0 m 2 =0 



which by eq (6.10) is also a sufficient stability condition 
for H( z i , z 2 ) . 

3 . Necessary Stability Conditions in Two-Dimensions 
By triangle inequality 




( 6 . 12 ) 



which, by eq (6.7) is equal to 
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mi=0 m 2 =0 





The right hand side of eq (6.13) converges if and only if 



Hi H 2 



1 1 

mi=0 m 2 = 0 



a 



mim 2 



< i 



(6.14) 



Equation (6.14) represents necessary conditions for stability 
of H( z i , z 2 ) . 

For the following sections and in Chapter VII, a 
mathematical generalization of the above conditions to 
N-dimensions , improved necessary conditions and a method to 
determine necessary and sufficient stability conditions in 
N-dimensions are presented using the theory of modern 
algebra as in Ref [53]. Notice from the unit sample response 
corresponding to the third order uncoupled transfer function, 
that the signs of a 0 i or ai 0 have no influence on the 
convergency behavior of <A> . This observation can be 
mathematically formulated in general as "symmetry condition" 
in the N-variable coefficient space. 

B. COEFFICIENT SPACE 

The expression (3j , . . . ,3^) denotes, using the notation 
established in Chapter IV, a fixed k-tuple of vector 



142 



indices, each of dimension N, i.e., for the third order 
uncoupled two-dimensional transfer function: 



(Si,5’ 2 ,S 9 ) = {(1,0); (0,1); (2,0)} . 

The characteristic equation Q( z”) , of an N-dimensional transfer 
function LCz) , where 



QC 



z) = 1 - L \ 



S. 

i 



(6.15) 



i=l 



is uniquely defined by the arbitrary k-tuple of real numbers: 
(Ai,...,Aj < .) and by the ordered k-tuple (S"i , . . . ,S^) of vector 
indices . 

For a fixed k-tuple of indices (31 , . . . ) the set of all 

k-tuples (Ai,...,A^) determines the coefficient space for the 

polynomial Q(z"). Each k-tuple or 'point' in the 

k 

k-dimensional Euclidian space IR is therefore associated 
with the polynomial of the form of eq (6.15). 

It is said that the point (Ai,...,A^) is a 'stable point' 
if QCz) associated with the k-tuple of real numbers corre- 
sponds to a stable transfer function. 

The set of all stable points in IR^ is called the 
'region of stability' in IR determined by the ordered 
k-tuple (Si , . . . ,3^) of vector indices. 



Example 6.1: 



The stability conditions for the second order coupled 
two-dimensional transfer function derived in the previous 
chapter define the region of stability in the three- 
dimensional coefficient space (A l5 A 2 ,A 3 ) as shown in 
Figs 6.1, 6.2, and 6.3. 

C. NECESSARY CONDITIONS 

In this paragraph, two necessary stability conditions 
in N-dimensions are derived, namely, that the absolute sum 
of the leading coefficients and of all coefficients in Q(z) 
must be strictly less than one. The term "leading 
coefficient" is defined in Appendix A. 

Theorem 6.1 : Let (A l5 ...,A^) be a point in the coefficient 

space and let A be the set of A^ which are leading coeffi- 
cients of the polynomial associated with (Aj , . . . ,A^) . If 
(Ajj.^jA^) is stable, then 



Proof: See Appendix M. 

Example 6.2 : 

The two-dimensional full quadratic transfer function 
has the form 




L(zi , z 2 ) 



1 



(6.16) 



2 



2 



1 — CtjoZl - &o i Z2 - OI11Z1Z2 



Ct20Zl 



CX02Z2 
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The set of leading coefficients is identified as 



A = {(1,1) ; (2,0) ; (0.2) } 

A necessary condition for stability of eq (6.16) requires 
that the absolute sum of the leading coefficients, which 
are an, a 20 , and a 0 a, is less than one, i.e., 

i cc 1 1 + a 20 + a 0 2 | < 1 • 

Another important result regarding necessary conditions for 
stability of an N-dimensional transfer function is: 

Theorem 6.2 : Let (ri,...,r^) be any N-tuple such that 

r^ = 0,1, or -1 for i = 1,...,N. Define the k-tuple 
(pi » • • • >p k > b y 



I> i = (r i , . . . , r N ) S x 

Then, if (Ai,...,A^) is a stable point, 

PiAi + ... + P k A k < 
Proof: See Appendix M. 



(6.17) 

it must be true that 
1 (6.18) 
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Example 6.3: 



For the full quadratic transfer function of the previous 
example, let (r l9 r 2 ) = (1,-1). Then the 5-tuple 
(Pi ,P 2 >P 3 »Pw >Ps ) is defined by 



P! = (1,-1) (1,0) 


= l 1 


(-1) 


p 2 = (1,-1) (0,1) 


= 1° 


(-1) 


_ _ ( -1 n s(l,l) 

P 3 - (1,-1) 


= -1 




= (i,-i) (2,0) 


= 1 




PS • (1,-1) (0 ’ 2) 


= 1 





Therefore, by Theorem 6 . 2 , if 

( A 1 , . . . , ) = (aio , . . . ,«02) 

is a stable point, then 

ot 1 0 - c*o 1 - 0.1 1 + CX 20 + ct 0 2 < 1 

Whenever all r- = 1 for i = 1,...,N the following corollary 
can be stated: 

Corollary 6.1 : If (Ai,...,A^) is a stable point, then 

Ai + ... + A k < 1 (6.19) 
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Theorem 6.3 : Let (ri,...,r N ) be any N-tuple such that 

= 0,1, or -1 for i = Define the k-tuple 

(pi > • • • 5 P k ) by 

S’. 

P^ = (ri,...,r N > 1 (6.20) 

Then, if (Ai,...,A^) is a stable point, it must be true that 
| p i A + ... + p ^A^ | < 1 (6.21) 

Proof: See Appendix M. 

Finally, following the above reasoning, corollary 6.2 is 
formulated as follows : 

Corollary 6.2 : If (Ai,...,A^) is a stable point, then 

| A i + ... + A^ | < 1 (6.22) 



D. SUFFICIENT CONDITIONS 

Theorem 6.4 : Let (Ai,...,A^) be a point such that 

k 

^ |Aj | < 1 (6.23) 

5=1 



then (Ai 




is stable. 



Proof: Follows directly by generalizing the derivation of 
section A2 to N-dimensions . 
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Example 6.4 : 

A sufficient condition for stability of L(z l3 z 2 ) as 
defined in eq (6.16) is: 

| ot 1 0 I + I OtQ 1 I + I ot 1 ! | + | Ot 2 0 | + I OC 0 2 | < 1 



A summarizing discussion of above stability conditions is 
presented in the next chapter. Note that necessary stability 
conditions and sufficient stability conditions are listed 
in tables 6.1 and 6.2 for two-dimensional filters up to 
sixth order and for selected N-dimensional filters, 
respectively . 

E. SYMMETRY CONDITIONS 

Theorem 6.5 : Let 5" = ( 6 1 ,..., 6^) be any complex vector such 

that 

i) |5^| = 1 for all i (6.24) 

ii) for i = l,...,k the_numbers p. defined by 

Si 

( 5 i > • • • > 5 
are all real. 

Then a point (A x , . . . jA^) is stable if and only if the point 
(AJ , . . . ,A£) is stable, where 

Aj = Pj^j f° r ^ = (6.26) 

Proof: See Appendix M. 



( 6 . 25 ) 

(6.25) 
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Example 6.4 : 

The third order uncoupled two-dimensional transfer 
function 

L(zi ,z 2 ) = 2 

1 — Otlo z l — &01 z 2 — 01 2 0 z 1 



If 6 = -1 



■ i = (-1,-1) (1,0) = -1 



P 2 



= (- 1 ,- 1 ) (0,1) = -1 



P 3 = (-1,-1) (2,0) = 1 



Thus, if (Aj ,A 2 5 A 3 ) is stable, then the necessary and 
sufficient condition of Theorem 6.4 prescribes that 
(-Ai,-A 2 ,A 3 ) is stable. This implies that the stability 
region of L(zi,z 2 ) is symmetric with respect to the (Aj) 
and (A 2 ) axis, which was already observed in Section A of 
this chapter. 



Transfer 

Function 


H(zi,z 2 ) = — — 

1 — aioZi “ cx o 1 Z 2 — 0(2 0 z 1 




STABILITY CONDITIONS 


Necessary 


| 0t2 0 | < 1 


Theorem 6.1 




|ai o| + |a 0 1 | + a 2 o <1 


Theorem 6.2 




| a 1 0 + oto 1 4 ot 2 0 1 < 1 
| — ot 1 0 + Oto 1 + 0i 2 o I <1 

|ai 0 - a 0 1 + a 20 1 <1 
|-aio - a 0 1 + a 20 | <1 


V Theorem 6 . 3 


Sufficient 


l ot 1 0 | + | Oto l| + I 0t 2 0 1 <1 


Theorem 6.H 



Table 6.1A: Sufficient and necessary stability conditions 

for the third order (uncoupled) two-dimensional 
transfer function. 
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H( z i , z 2 ) = — 

3 3 

, T 7 a -mi -m 2 
1 - L L mim 2 Zi a z 2 

mi=0 m 2 =0 

(0<mi+m 2 < 3) 



STABILITY CONDITIONS 



Necessary |a3o + <*21 + «i 2 + a 0 3| < 1 

|aio + ao 1 + Qtn + a 20 + cto2 + a 2 o + <^12 + a 3o + 0103 ! 


Theorem 

6.1 

\ 


< 1 




| — ct 1 0 + a 0 j + an - a 2 o + ao2 + a 2 i + a i2 - a 30 + a o3| 




< 1 

I Ot 1 0 - a 0 1 - a ll + a 20 + a 02 - a 21 + 0ti 2 + a30 - ao 3 1 


> Theorem 

6.3 


< 1 




|-aio - a 0 1 + an + a 20 + a 0 2 - ot 2 1 - ai 2 - a 3 o - ao3| 




< 1 


J 


3 3 




Sufficient I l la^J < 1 

mi=0 m 2 =0 

(0<mi+m 2 _< 3) 


Theorem 

6.4 



Transfer 

Function 



Table 6. IB: Sufficient and necessary stability conditions 

for sixth order (coupled) two-dimensional 
transfer function. 
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Transfer Function 

Hu , ,10 - 



^ z. z 



i- _ Z1 H. H_ ro z wrij ' O 2$ 

ft, ,'o M,*= o 

to * n, + nr. 1+ m i 61) 



STABILITY CONDITIONS 



Necessary : 

' ^ + <*.,.+ «,» + ^ <+ <> U1 
^ ^lt>9 + ^\a °fj90 ^Zoj * i<5»9 + ^iu ~ °M0| * ^ C{1 ^ > 

^ ^IM ~ Sip + <Af, ^ +■ «< 0?o + <S ,1 " <^'o T ^ 1 ®, - ^ i ^ * 

f 

^ °Sj +■ ~ 4* ^1*9 4 “H Hjtx 4- ~ °/lc, - S r> I ^ i 

^ 0,1:5 S, ^ “H + ^0 L 4“ ^3 “ 'W ~ ^ ^ 

I ^ + <4®„ - c/ /?p + 4- <\ u V ci onl - 4 *»v ~ <**¥ )^ 1 

I C **iD) "* <4l3 ' ‘SpO ** S*, 4- ^<310+* ^ - q, ’ ,, ^ ^ “ 



Sufficient 

■/- 7 



['[ [ 

tt^'O Hr .^0 VHjrO 

6° 4 ^1+ rv,^ nij tv> 



Theorem 6.1 



Theorem 6 . 3 



Theorem 6.4 



Table 6.2: Sufficient and necessary stability conditions 

for full quadratic three dimensional transfer 
function . 
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FIG 6.1: STABILITY REGION OF BILINEAR FILTER 
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FIG 6.2: STABILITY REGION OF BILINEAR FILTER 
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FIG. 6.3: STABILITY REGION OF BILINEAR FILTER 




VII. NECESSARY AND SUFFICIENT CONDITIONS 



IN N- DIMENSIONS 



A. INTRODUCTION 

A stability region is said to be convex if the loci of 
all points connecting two stable points lies inside the 
stable region. It is well known that the stability area 
for a second order one-dimensional transfer function satisfies 
convexity, but it can be 'shown that the stability region of 
the third order one-dimensional transfer function is not 
convex. It is, therefore, conjectured that the bounding 
surfaces of a multi-dimensional higher order filter stability 
region cannot be described in terms of plane surfaces which 
explains why the following method to derive necessary and 
sufficient stability conditions is rather complicated to 
apply. The advantage of this approach lies in the fact that 
the method must be applied only once for a particular type 
of transfer function, which results in a set of stability 
conditions. These conditions can be easily applied by the 
design engineer. 

B. METHOD TO DERIVE NECESSARY AND 

SUFFICIENT STABILITY CONDITIONS 

The method to derive necessary and sufficient stability 
conditions consists of two steps: 

1. Derive necessary and sufficient conditions for one 
variable polynomial with complex coefficients using 
Rouche 1 s theorem . 
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2. Apply the following theorem and corollary to state 

necessary and sufficient conditions for N-dimensions . 



Theorem 7 : Let (Si,..., 6^) be any N-tuple of complex 

numbers such that r = max |S^| <_ 1 . Then Q(zi,...,z^) 
is stable if and only if for all 6 the polynomial 



q^(z) = Q ( z rp 




(7.1) 



is stable. 

Proof: See Appendix N. 

Corollary 7 : The polynomial Q(z") corresponds to a stable 

transfer function if and only if for all complex numbers 3T 
such that 



max |S i | = 1 , q^(z) = Q( zS i , . . . ,zS N ) 



(7.2) 



is stable. 



Proof: See Appendix N. 



C. EXAMPLE 7.1: FIRST ORDER TRANSFER FUNCTION 

IN N-DIMENSIONS 

Using the notation of the previous chapter the first 
order case in N-dimensions is defined as follows : 



Q(z) 



N 




i=l 



(7.3) 
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where 



Ai 

A2 



(1,0,. . . ,0) 
■(0,1,. . . ,0) 



N 



( 0 , 0 ,. .. , 1 ) 



The one variable polynomial 



N 



q(z) = 1 - ) A^6^z for max |6^| = 1 



i=l 



is known to be stable if and only if 



N 

7 A. 6 . 

Z_ 11 

i=i 



< 1 



Thus by corollary 7 Q(z) is stable if and only i 
holds for all 6 , where max 1 6 - I = 1 . For 6. 
eq (7.4) results in 

|Ai| + ... + | A n I < 1 



which is a necessary condition. But by Theorem 6 



(7.4) 

f eq (7.4) 

= ±1 

(7.5) 

4 , eq (7.5) 
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is also a sufficient condition. In summary the first order 
N-dimensional filter is stable if and only if the sum of 
the absolute coefficients is less than one. 

D. EXAMPLE 7.2: N-DIMENSIONAL TRANSFER FUNCTION WITH 

POSITIVE COEFFICIENT DENOMINATOR POLYNOMIAL 

If Q(z") is defined as: 
k 




i=l 



and all a^- >_ 0 , then a necessary and sufficient condition 

i 

for stability is that the sum of all coefficients is less 
than one. 

This follows immediately by applying Theorem 6.2 for 
r^ = 1 to derive the necessary condition and Theorem 6.4 
for the sufficient conditions, which for all coefficients 
positive are identical. Note that the stability condition 
for the above cases as well as for the bilinear case 
(without proof) are tabulated in the following Table. (Table 7). 

E. SUMMARY 

Several powerful methods have been developed to solve 
the stability question of N-dimensional recursive digital 
filter. Ultraspherical polynomial and derivative operator 
methods were used to derive stability conditions for specific 
low-order two-dimensional cases. Methods for the derivation 
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of sufficient and necessary stability conditions in two- 
dimensions and in N-dimensions and for necessary and suffi- 
cient stability conditions in N-dimensions were presented 
and applied to compute low order two-, and N-dimensional 
stability conditions. 

A more complete derivation of stability conditions for 
low order two- and N-dimensional recursive digital filter 
and an arrangement in tabular form are left for future 
research. With the methods presented to solve the stability 
question the second important part, i.e., approximation 
of a given specification by a transfer function, will be 
investigated in the following chapters. 
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VIII. PROPAGATION OF TRANSFER FUNCTION COEFFICIENTS IN 



THE UNIT SAMPLE RESPONSE OF N-DIMENSIONAL FILTER 
A. INTRODUCTION 

In the previous chapters, stability criteria for 
N-dimensional filters and implementation forms were developed 
which will be applied in the following chapters to derive 
time domain design techniques of N-dimensional recursive 
digital filter. The design of a recursive digital 
filter is achieved in two steps, (1) choose filter coeffi- 
cients to approximate a given response, and (2) ensure that 
the resulting filter is stable. In Chapters IV, V, VI, and 
VII, several techniques to solve the stability question were 
developed and stability criteria for N-dimensional and 
specific two-dimensional recursive filter derived. In the 
following chapter the nonrecursive combinatorial formula of 
Theorem 4.4 will be utilized to relate transfer function 
coefficients with the entries of the unit sample response 
in N-dimensions . It will be specifically shown that the 
coefficient of the characteristic equation a^— ^ appears 
only once in linear form and with coefficient one in the 
coefficient h(n) of the corresponding unit sample response. 

It also appears in nonlinear form, with non-unity coeffi- 
cients in all h(j), where n < J . 

This important correspondence will be used in Chapter IX 
to construct an algorithm for the extraction of an all-pole 
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transfer function from a given unit sample response. The 
extraction of a rational N-dimensional transfer function will 
be achieved in the same chapter utilizing a method which is 

r 

essentially a Pade approximant extended to N-variables . 

Both of these methods will be incorporated in an algorithm 
extracting a parallel arrangement of low order sections. 

In Chapter X the time domain design step 2 will be 
achieved by approximating a given unit sample response 
(h 0 (n)} , where 0 _< n.^ 00 » 0 < i _< N by {h(n)} , 

where 0 _< <_ ; and then using transfer function 

extraction methods. In step 2 of the design procedure, the 
approximating transfer function is tested for stability. The 
quality of the time domain design is evaluated comparing 
the spectrum of {h(n)} to the spectrum of the extracted 
transfer function. 

B. PROPAGATION RULES IN N-DIMENSIONS 

It has been stated in Theorem 4.5 that the Taylor 
expansion of 



H(z) = , 

1 - Q( z ~ ± ) 



where 




x 

( 0<i i <°°) 
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has the form: 



HCz) 




— x -n 



h( n) z 



( 0 <ri£<°°) 



The coefficients h(n) are related explicitly to the a ('j') 
by theorems 4.4 and 4.5 as follows: 



h(n) = 




■to > 
b o 






( 8 . 1 ) 



where 






( 8 . 2 ) 



In the next paragraphs , three cases for the appearance of 

a particular coefficient a r ^- x , where (0 < p < k) , 

P 

in the coefficient expression h(n) of the Taylor expansion 

of H(z) are investigated, namely (S’ ) > (n) , (S’ ) = (n) , 

¥ ¥ 

and (S ) < (n) . 

P 
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1. Propagation Rule for a,*- . , where (S ) > (n) 

p 



The notation (S' ) > (n) is defined as 

P 



(S_ 


) 


> (n ) 


p 1 


l 


(S 


) 


> (n ) 


P 2 


2 



or 



(S ) 
p N 



> <V 



For example, for (S ,S p2 ) = (2,1) the inequality 
(2,1) > (ni,n 2 ) designates the shaded area in the 
following graph: 
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of the N-variable 



Propagation Rule I : The coefficient a^g- ^ 

characteristic equation of H("z) does not appear in h(n) , 



whenever (S ) > (n) , where h(n) is the n 

P 

of the’ Taylor expansion of H(~z) . 



-th 



coefficient 



Example 8.1: 
For 



H( Zi , z, ) = 

1 z -1 -1 -3 - 5 “ ** 

1 - (X 2 0 Z 1 — Ct 0 1 z 2 — 01 3 0 z 1 — 0( 5 4 Z j Z ? 



( 8 . 6 ) 



the corresponding unit sample response {h(nj,n 2 )} is 
computed from eq (4.21). In the case where (n l 5 n 2 ) < (5,4) 
eq ( 8 . 2 ), which becomes for the above transfer function 

t (l,0)‘ (l5 ° ) + t (0,l)‘ (0 ' 1) + t (3,0)' (3,0) + t (5,4)* (5,4) = (ni,n2) 
can only be solved, if and only if, ^(5 4 ) = ^ • 



2. Propagation Rules for a^g ^ , where (S^) = ^ 



The notation 


“I 

II 


(n) 


is 




( s _ ) 
Pi 


- 


(n x ) 


and 


(S ) 
p 2 

• 


= 


(n 2 ) 


and 


• 

• 

(s^ > 
Pn 


= 


( V 
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If (Sp) = (n) , then there exists one term in h(n) where 

t^- equals one and all other t's are identically zero, 

' P 

therefore eq (8.1) can be rewritten as follows: 



h(n) = 



v° 

ts- , 0,...,0 
s p 



ts. 



a (S ) 



terms of 
order < S_ 



t S_ =1 



, , terms of 

a (S ) 

p ; \ order < S_ 



In this case a,s- . appears in linear form with coefficient 

C P 

one. This observation is stated as follows: 

Propagation Rule II : The coefficient a,g- ^ of the N-variable 

characteristic equation of H(z) , appears in h(n) in linear 
form with coefficient one, whenever (Sp) = • 

Example 8.2 : 

The coefficient a^ ^ = a 54 of the same transfer 
function as in example 8.1 appears in h(5,4) in linear 
form with coefficient one, since 



t (l,0) (1,0) + t (0,l) (0,1) + t (3,0) (3,0) + t (5,4) (5,4) = (5,4) 

whenever t^g ^ = 1 and all other t’s are zero. The terms 
of order < (5,4) are identified as: 
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^ 1 , 0 ) = 2 
and 

t (l»0) = 5 



t (0,l) " 4 ’ t (3,0) = 1 ; t (5,4) 

t (0,l) = 4 ’ t (3,0) = 0 > 



Therefore, using eq (4.48a) 

h(5,4) = as*t + 105 aio ad a 3 o + 112 ad aoi 

If the I^ 1 diagonal cut of the unit sample response is 
defined to be the set of indices (n) = (m,...,n^) , such 

that 




(8.7) 



then the following corollary to Propagation Rule II can be 
stated : 

Propagation Rule IIA : Consider the first diagonal cut of 

{h(n)} , then each entry corresponds directly to a first 

order coefficient a,«- \ of the N-dimensional characteristic 

( P 

equation of H(z) . In fact, for 
N 

[ “i = 1 ’ 

i=l 



h(n) 




whenever (S ) = (n) 
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Example 8.3 : 

The unit sample response entries h(n) along the first 
diagonal cut corresponding to the transfer function 



H(zi,z 2 ) = 



1 -CtiooZl - aoioZz 1 - CtooiZs 1 “ C t 110Zi 1 Z2 1 ~ as 4 3Z1 5 Z 2 **Z3 3 



are h ( 1 , 0 , 0 ) = aioo 

h(0 ,1,0) = ao 10 
h( 0 , 0 , 1) = ao 0 1 



This follows directly from eq (8.2), which can be written 
as follows 



t (l,0,0) (1,0,0) + t (0,l,0) (0,1,0) + 



+ t (l,l,0) (1,1,0) + t (5,4,3) (5,4 ' 3) 



= {( 1 , 0 , 0 ); ( 0 , 1 , 0 ) ; ( 0 , 0 , 1 ) } . 



(End of example) 



W 

It is also observed from eq (8.2), that a,g- . appears 

. . _ K 

with coefficient one, whenever (n) = W(S ) , then 

P 
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h(n) = 




% =W 
P 



( terms of 

order < W (S’ ) 
P 



= a 



W 

( V 




Therefore, the following rule can be stated: 

Propagation Rule III : The coefficient a . of the 

_ P ; 

N- variable characteristic equation of H(z) appears in h(n) 

"t h _ 

to the W L power with coefficient one, whenever (n) = W(S^) . 

Example 8.4 : 

The coefficient a(5,4) = a 5 i* of previously used 
transfer function appear in h(10,8) with coefficient one 
to the second power, which can be identified from 



( 1 , 0 ) 



+ t 



( 0 , 1 ) 



(0,1) + t (3 0) (3,0) + t 



(5,4) 



(5,4) = (10,8) 



where for t^g ^ = 2 all t's must be zero. The terms of 
order < (10,8) are computed from: 
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t (l: 


,0) 


= 10 


t (0,l) 


= 8 


t ( 1 : 


,0) 


= 7 


t (0,l) 


= 8 


t (l ; 


,0) 


= 4 


t (0 ,1) 


= 8 


t (l ; 


,0) 


= 1 


t (0,l) 


= 8 


t (l : 


,0) 


= 2 


t (o ,1) 


= 4 


t (l ; 


,0) 


= 5 


t (0 ,1) 


= 4 



Thus 



2 18 T 10 s 
h( 10 , 8 ) = asi. + jgT'^T a ioa 0 1 



14 T 4 8 2 

g | ^ i 2 1 ^10^X01^X30 



8! 2 4 

2] "4; aioao 10x30054 



t ( 3,0) 


= 0 


t (5,4) 


t (3,0) 


= 1 


t (5,4) 


t (3,0) 


= 2 


t (5,4) 


t ( 3,0) 


= 3 


t ( 5 , 4 ) 


t (3,0) 


= 1 


t ( 5 , 4 ) = 


t (3,0) 


= 0 


t (5,4) 


16! 


7 8 


1 



~ 7 ! ~' 8! ~ a io a oi a3 o 



12! 8 3 

+ gT - jT a ioaoia2o 

10! 5 4 

g i ^ y aiooto 1CX54 



"th _ 

Propagation Rule IV : In the I diagonal cut of h(n) , 

"th 

each entry contains an I L order coefficient a,-^ <. of the 

^ P 

N-dimensional characteristic equation of H(z). In fact, 
h(n) = a^g- ^ + (terms of order < S^) , where (S ) = (n) . 
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where (S ) < (n) 
P 



as 



3. Propagation Rules for a,^- \ 5 

c v 



The notation (S ) < (n) is defined (as in Theorem 4.1) 

P 



(S_ ) < (m) , 
P 1 — 

and (S ) < (n 2 ) 

P 2 “ 



and (Sp ) < (n^) , where the case 

(Sp) = (n) is explicitly excluded. 



For example, for (S jSp^) = (2,1) the inequality 

(2,1) < (r.i ,n 2 ) designates the area, which is shaded in the 

following graph: 



n 2 = 0 1 2 3 4 
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Example 8.5 : 

The set of coefficients {a.. . , } , where ii + i 2 = 2 , 

Ui ,i 2 ) 

appear in the coefficients of the Taylor expansion of 



H(Z!,Z 2 ) = 

-1 -1 -1 -1 -2 -2 -5 -■» 

1 “ C6lO z l * CXq 1 z 2 “ anZi Z 2 “ a 2 0 Z l ” CXo 2 z 2 ” CX54Z1 z 2 



( 8 . 8 ) 



along the second diagonal cut, i.e., in {h(n 15 n 2 )} , such 
that ni + n 2 = 2 , since 

t (l,0) (1,0) +t (0,l) (0,1) +t (l,l) (1,1) +t (2,0) (2,0) +t (0,2) (0,2) 

+ t (5,4) (5 ’ 4) = (2 ’ 0) 

t (l,0) (1,0) + t (0,l) (0,1) + t (l,l) (1,1) +t (2,0) (2,0) +t (0,2) (0,2) 

+ t (5,4) (5,4) = (1,1) 

t (l,0) (1,0) +t (0,l) ( ° 5l) +t (l,l) (1,1) +t (2,0) (2,0) +t (0,2) (0,2) 

+ 1(5,4) (5,4) = 

Consequently, we obtain by inspection 



h( 2 , 0 ) = a 2 0 


+ 


2 

<*i 0 


h( 1 , 1) = ail 


+ 


2c*2 0 do 1 


h( 0 , 2 ) = ao 2 


+ 


2 

ao 1 
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For (S^) < (n) , the following rules can be stated: 



Propagation Rule V : The coefficient a^g- ^ of the N-variable 

characteristic equation of H(z”) appears in h(n) for 
('Sp) < (n) , whenever equation (8.2) is satisfied. 

W 

Propagation Rule VI : The coefficient a^ ^ of the 

P 

N-variable characteristic equation of H(z) appears in h(n) 
for W(S" ) < (n) , whenever equation (8.2) is satisfied. 

p 



C. PROPAGATION DIAGRAMS FOR TWO-DIMENSIONAL TRANSFER 
FUNCTIONS UP TO THE SIXTH ORDER COUPLED CASE 

The Taylor expansion coefficients of two-dimensional 

sixth order coupled case, i.e., 



’ _i - 1 -1-1 -2 -2 -3 

1 - Ot i o Z i - CX o 1^2 “ CX i iZj Z 2 “ (X 20 Z 1 " CXo 2 Z 2 0^30^1 

( 8 . 9 ) 

— -2 -1 -1 -2 II 

- a 2 izi z 2 - ai 2 zi z 2 - a 0 3 Z 2 

will be investigated in this section. It was shown in Section 
D of Chapter IV, that in two-dimensions the unit sample 
response matrix A is obtained by writing the Taylor 
expansion of H(zi,z 2 ) as matrix product in the following 
manner 



H(zi,z 2 ) = (1 z i 1 z7 2 zT 3 ...) A (1 z 2 1 z 2 2 z 2 3 •'•) t 



17 ^ 



where 



1 


h(0 ,1) 


h( 0 , 2 ) 


h(0 ,3) 


h( 1 , 0 ) 


h( 1 , 1) 


h( 1 , 2 ) 


h( 1 , 3) 


O 

n 

CM 

V-/ 

A 


h( 2 , 1) 


h(2 ,2) 


h( 2 , 3 ) 


h( 3,0) 


h( 3,1) 


h( 3 , 2 ) 


h( 3,3) 



If for a moment the case is considered, where in eq (8.9) 
all coefficients are zero, except aio and ao i , i.e., 



H( z ! , z 2 ) = ^ — (8.10) 

1 ~ aiozj ~ 0ioi z 2 

then the corresponding unit sample response matrix A is 
shown in Fig (8.1). The coefficients aio and aoi > which 
are by definition (see Appendix A) also the leading coeffi- 
cients of eq(8.10), appear by propagation rule IIA in h(l,0) 
and h(0,l) in linear form with coefficient one, and 
"propagate" by propagation rule V for all (n^n^) > (1,0) 
and (nj,n 2 ) > (0,1) , respectively. These considerations 

are easily identified in Fig 8.1. 

An abstract method to graph A is shown in Fig 8. , 

in which each entry of the unit sample response matrix is 
replaced by a cross, and single lines outline the sector of 
A in which the coefficients cti 0 and a 0 i "propagate", 
i.e., appear whenever eq (8.2) is satisfied. 
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Using similar symbology the leading coefficient 
propagation of the third-order uncoupled transfer function, 

i.e. , 



H( z i , z 2 ) 



1 

- 1 - 1 - 2 
1 — Cti o Z i - OI 01 Z 2 “ 0t20 z l 



( 8 . 11 ) 



. W 

is shown in Fig 8.2a. The propagation of a 20 , which 

follows from propagation rule VI is outlined and identified 

by the coefficient CX 20 raised to the W power. 

Figure 8.2c shows the leading coefficient propagation 

for the second-order coupled (or bilinear) transfer function. 

W 

It is observed that aj 1 propagates m accordance with 
propagation rule VI from the w"*^ 1 entry along the main 
diagonal for all (ni,n 2 ) > (1,1). 

Figures 8. 2d and 8.3e display the leading coefficient 
propagation of the fourth-order coupled (or full quadratic) 
and the sixth order coupled transfer function, which follows 
from propagation rule IV, where 1=2, 3, respectively, 
and propagation rule V. 

The following table relates two-dimensional transfer 
function and leading coefficient propagation charts up to 
the sixth order coupled case. 

D. SUMMARY 

The six propagation rules developed in Section B from 
the non-recursive combinatorial formula of Theorem 4.4 
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establish a relationship between the coefficients of a 
transfer function in N-dimensions to the coefficients of 
an N-dimensional unit sample response. The rules state 
explicitly that a transfer function coefficient a (J) 
appears only in h(T) in linear form and with coefficient 
one and propagates in all h(n) , where (n) > (i) . 

Conversely , if from a given unit sample response it is 
desired to identify the transfer function coefficient > 

then by propagation rule II entry h(i) contains a ('J') > 

as well as nonlinear terms composed of a (J) > such that 
j < i . But, if the unit sample response of the transfer 
function H^(zi,z 2 ) having {a^j^} coefficients, where 
j < T , is known then a (J") is equal to the difference 
of h^^ and the i "*"^ 1 entry of the unit sample response 
corresponding to H^(zi,z 2 ). Moreover, it will be shown 
in part B of the next chapter that this idea can be imple- 
mented to extract all a (T) coefficients, where 
N 



I 

i=l 



at the same time. 
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TABLE - 9 ; EXAMPLES OF LEADING COEFFICIENT PROPAGATION IN TWO DIMENSIONS 
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otoi do i do 1 
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FIG 8.1a UNIT PULSE RESPONSE A OF SECOND ORDER UNCOUPLED FILTER 
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FIG 8.1b.: PROPAGATION OF a 10 AND a 0 i COEFFICIENTS 
(SECOND ORDER UNCOUPLED FILTER) 
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FIG 8.2c-. PROPAGATION OF On COEFFICIENT 
(SECOND ORDER COUPLED FILTER) 
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.2d: PROPAGATION OF a 20 ,ctii,ao 2 COEFFICIENTS 
(FOURTH ORDER COUPLED FILTER) 
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(FIFTH ORDER COUPLED FILTER) 
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FIG 8.3d: PROPAGATION OF a 0 3 COEFFICIENT 
(FIFTH ORDER COUPLED FILTER) 
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(SIXTH ORDER COUPLED FILTER) 



IX. EXTRACTION OF TRANSFER FUNCTIONS FROM A 



GIVEN UNIT SAMPLE RESPONSE 
A. INTRODUCTION 

It was outlined in Section D of the previous chapter 
how the propagation rules can be used to identify a transfer 
function coefficient from a given unit sample response. 

This technique is elucidated and is used to develop an 
algorithm which extracts an all-pole transfer function from 
a unit sample response in N dimensions. 

In Section B, an extraction method for rational transfer 
functions is presented which is essentially a Pade approxi- 
mate generalized to N-dimensions . This procedure is shown 
to require the simultaneous solution of a number of linear 
difference equations, which will preferably be done using a 
general purpose digital computer. 

Finally, in Section E, an approximate partial fraction 
expansion in N-variables in terms of all-pole and/or rational 
low order transfer functions is presented. 

A measure of quality for the extraction methods will be 
established in the next chapter, in which the above tech- 
niques are incorporated to solve step 1 of the time domain 
design procedure. 
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B. EXTRACTION OF AN N- DIMENSIONAL ALL-POLE TRANSFER 
FUNCTION FROM A GIVEN UNIT SAMPLE RESPONSE 

The given N-dimensional unit sample response (h(n) } , 

where 0 _< n^ £ 00 and 0 < i _< N is chosen to be normalized, 

i.e., h(0) = 1 . The transfer function coefficients having 

first order, i.e., the set of a,T. , for which the sum of 

all indices of (T) is identically equal to one, for example, 



[a (l,0,...,0) ;a (0,l,0,...,0)’‘ ,, ’ a t0,...,0,l,0)’ a (0,...,0,l) 1 

are known by propagation rule II A to be identical to the 
corresponding set of unit sample response entries h(n) , 
taken along the first diagonal cut. 

The entries {h(n)} , where (n) is taken along the 
second diagonal cut contains, by propagation rule IV, all 
second order coefficients, in linear form with unit coeffi- 
cients, as well as first order nonlinear terms. To identify 
the set of second order coefficients, i.e., ^- a ('f)^ 2 ’ 

the unit sample response entries {hi (n) } corresponding to 
the transfer function Hi(z) , which has only first order 
coefficients, is computed and the subtraction 
{Ah 2 (n)} = (h(n)} - {hiCn)} is performed for every term 
along the second diagonal cut. The resulting set of 
coefficients Aha(n) is identical to the corresponding 
{ a C i ) } 2 • It is apparent that this technique leads to an 
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algorithmic approach in which the order coefficients, 

i.e., ^ a (i)^ are identified from the set of differences 



Ah^(n) = h(n) - h^_^(n) 



taken along the k"^ diagonal cut. If this procedure is 
continued until (Ahj(n)} = {(0)} , then the all-pole 
transformation Cz) is the exact transfer function 
corresponding to the given unit sample respone (h(n)} . 

It will be shown in the next section that the above 
technique can be quite readily implemented on a general 
purpose or table top computer. It is noted that the order 
of the extracted transfer function in each variable is less 
than or equal to the size of the unit sample response along 
the corresponding mesh axis. 



C. COMPUTER ALGORITHM FOR THE EXTRACTION OF AN 
N- DIMENSIONAL ALL-POLE TRANSFER FUNCTION FROM 
A GIVEN UNIT SAMPLE RESPONSE 

The following algorithm incorporates, in addition to the 
concept developed in Section B, a mean square error 
criterion which will be used to terminate computations. 

Example 9.1 ; Extraction of a two-dimensional digital 
transfer function from a given unit 
pulse response 

The N-dimensional algorithm, which is outlined in block 
diagram form in Fig 9.1, has been implemented on an HP 9830A. 
The next pages contain an example which is designed to 
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FIG 9.1 EXTRACTION OF AN N-DIMENSIONAL 
ALL-POLE TRANSFER FUNCTION 
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illustrate the application of the all-pole transfer function 
extraction method to a two-dimensional unit pulse response 
which is chosen for convenience to be normalized, i.e., 
h(0,0) = 1 . The unit pulse difference matrix and the 
mean-squared-error are listed for every traversion of the 
main loop. 



Loop 1 ♦ Step 1: (Notice that the step number corresponds 

to the block number in Fig 9.1). The given unit pulse 
response {h(ni ,n 2 ) } is listed in Table 9.1, where 
0 ni <_ 10 , 0 <_ n 2 £ 10 . The algorithm will be 

_ 7 

terminated if the mean squared error is less than 10 , 

_ n 

i.e. £ = 10 . A graphical representation of {h(m ,n 2 ) 1 

is shown in Fig 9.2. Step 2: Initialize D(zi,Z 2 ) , I , 

{a^— , and (Ahi(ni,n 2 )} = (h(ni,n 2 )} . Steps 3,4: The 

coefficients ^ and a^ ^ are extracted along the 

first diagonal, i.e. nj+n 2 = 1 of {Ah^n! ,n 2 ) } , where 



a (l,0) 



. 52 



a (o,i) 



.65 



and 



The two-dimensional characteristic equation is identified 
as D(zj,z 2 ) = 1 - .52 zT 1 - .65 z 2 * . Step 5: The unit 

pulse response {h 1 (n 1 ,n 2 )} corresponding to 
H 1 (z 1 ,z 2 ) = D ■''(z 1 ,z 2 ) is computed using the recursive 
technique of theorem 4.1. Step 6: The difference matrix 

{AhjCn! ,n 2 ) } is obtained by matrix subtraction: 



19 ^ 






{Ah 2 (n ls n 2 )} = (h(n 1 } n 2 ) - {h 1 (n 1 ,n 2 )> and is listed in 

table 9.2. The mean square error e is computed as follows: 

10 10 
£ . V 

nj=0 n 2 =0 

_ 7 

Steps 7,8,9: The tests e versus 10 and I versus 20 are 

traversed and the main loop index I incremented to 2. 

Loop 2 . Steps 3,4: The coefficients 

a (2 0 ) = Ah 2 < 2 , 0 ) = -.11 
a (l 1) = Ah 2 ( 1 ,1) = “ • 2 
a (0 2 ) = Ah 2 ( 0 , 2 ) = .003 



I 



[Ah 2 (n, ,n 2 ) ] 



29. 81 



are identified off the second diagonal cut, i.e., n x +n 2 = 2 , 
of {Ah 2 (n 15 n 2 )} . The two-dimensional characteristic 
equation becomes 

D(zi,z 2 ) = 1 - .52Z7 1 - .65z 2 J + .llz7 2 + .2z7*z 2 1 - .003z 2 2 . 

Steps 5,6: {h 2 (nj,n 2 )} is calculated and {Ah 3 (nj,n 2 )} . 

(see table 9.3), is formed. e = 9.854E -3 . Step 7 to 9 : 

Tests traversed and I incremented to I = 3 . 

Loop 3 . Steps 3,4: The coefficients a^ ^ ^ = Ah 3 (n ls n 2 ) 

for n 3 +n 2 = 3 are identified as 
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a ( 3 , 0 ) 


= Ah 3 (3,0) 


= .034 


a ( 2 , 1) 


= Ah 3 ( 2 , 1) 


= .0032 


a ( 1 , 2 ) 


a ( 0 , 3 ) 


= 0. 



The 
D( Zj ) 



two-dimensional 



z 2 ) = 1 - 0.52Z! 



-l 



- 0.034z 



-3 

1 



characteristic 

- 0.65Z” 1 + .llzj 2 
- 0.0032 z 1 2 z 2 1 



equation becomes 
+ .2z - 0.003 z 2 2 



Steps 7,8,0: Tests traversed, 1=4. 

Loop 4 . Steps 3,4: The coefficients are read from the 

fourth diagonal of {Ah 4 (n x ,n 2 ) } as 



a ( 4 , 0 ) * 027 

a ( 1 , 3 ) = " ,053 



a ( 3,1) ’ a ( 2 , 2 ) = a ( 0 , 4 ) ” 0 



The two-dimensional characteristic equation becomes 



D(zj ,z 2 ) 



1 — 0 . 52zj — 0.65 z 2 t . llz i *2z j z 2 — 0.003z 2 

- 0.034z7 3 - 0.0032z7 2 z 2 1 - 0.027z7 4 + O.OSSz^^i 3 



Steps 5,6: The unit pulse difference matrix (Ah 5 (ni,n 2 )} 

is listed in table 9.5. e = 5.07566E -3 . 

Steps 7,8,0: Tests; 1=5. 
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Loop 5 » 



Steps 3,4: The coefficients along diagonal five 



are : 

a (4 1) = 0.0351 

a ( 2 , 3) = -0-0107 

a ( 5 , 0 ) " a (3,2) = a C 1 , 4 ) = a ( 0 , 5) = 0 ’ 

and 

D(z ,z ) = 1 - 0.52 zi ! - 0.65Z2 1 + .llzT 2 + .2zi 3 Z2 2 - .003Z2 1 

- 0.034zi 3 - 0.0032 zi 2 Z2 1 - .027zi 4 + •053 zi 1 Z2 3 

- 0.0351zT‘*z7 1 + 0.0107zT 2 zI 3 

Steps 5,6: {Ah 6 (ni,n 2 )} is listed in table 8.6. 

e = 1.71*10 ^ . Steps 7,8,9: Tests; 1=6. 

Loop 6 : Steps 3,4: The coefficient a^ g ^ = -0.067 is 

the only non-zero entry along the sixth diagonal of 
(Ah 6 (n 1 ,n 2 )} . The two-dimensional characteristic equation 
becomes : 

D( Zi , Z 2 ) = 1 - 0.52zj - 0.65Z2 + .llzi + .2zi Z 2 - 0.003Z2 

- 0.034zi - 0.0032 zi 2 Z2 1 - 0.027zi 4 + 0.053z~ 1 Z2 3 

- 0.0351zT 4 zI 1 + .0107zT 2 zI 3 + 0.067zT 6 
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